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BUCKLING OF CONTINUOUS BEAMS ON ELASTIC 
SUPPORTS 


By J. M. KLITCHIEFF (High School of Technology, Belgrade) 


[Received 3 January 1949] 


SUMMARY 

It is sometimes necessary to select the rigidity of intermediate elastic supports 
of a beam so that they behave as though they were absolutely rigid when the beam 
buckles on being subjected to longitudinal compression. The usual method of solving 
the problem leads to the tedious work of developing a determinant of order m— 1 and 
of solving an equation of degree m—1, m being the number of spans. The author 
represents the deflexion curve by trigonometric series and arrives at equation (6’), 
which is more convenient for practical use. The value of the rigidity obtained differs 
by less than 2-6 per cent. from the true value. 


Introduction 

Let us consider a continuous beam (Fig. 1) of constant cross-section 
simply supported at the ends on rigid supports 0 and m and having equally 
spaced intermediate elastic supports 1, 2,..., m—1 of equal rigidity. The 
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beam is subjected to longitudinal compression by forces P applied at the 
ends. It is sometimes necessary te select the common rigidity of the 
intermediate supports so that they will not deflect when the beam buckles 
and therefore be equivalent to absolutely rigid supports. 

This problem was first discussed by I. G. Boobnov (1), who used the 
following method. Calculating the deflexion of a simple beam submitted 
to the simultaneous action of compressive force P and lateral reactions 
B,, R,,..., R,,-, of the elastic supports, and putting the deflexion of a 
support 7 equal to UR, (W being the ‘spring constant’ of the support, that 
is, the deflexion produced in the support by the unit load), he arrived at 
m—1 equations with 7 a: —l, 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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The buckling of the beam becomes possible when these linear and 
homogeneous equations yield a solution for R,, R,,..., R,,, different from 
zero. Hence the critical values of P are found by setting the determinant 
of the equations equal to zero. It leads to an equation of degree m—|] 
and the roots determine the values of P corresponding to the buckling 
of the beam into 1, 2...., m—1 half-waves (shown on Figs. 2 and 3 for 
m = 3). If Wis a given number some of these values can become smaller 





than the force P,, corresponding to the buckling into m half-waves (Fig. 4, 
for m = 3), when every intermediate support becomes an inflexion point. 
To determine the limiting value of % at which the intermediate supports 
still behave as though they were absolutely rigid, I. G. Boobnov introduced 
into the equation for P the value P,, and arrived at m—1 values of Y as 
roots of the equation. The smallest of them is the limiting value required. 
This method of solving the problem is a logical one, but it leads to 
developing a determinant of order m—1 and to solving an equation of 
degree m—1. This work becomes tedious if the number of spans is not 
small. It was done by I. G. Boobnov for m = 2, 3, 4, 5, 6, 7, 9, and 11. 
From the numerical values obtained he concluded that % approaches 
a certain limit as the number of spans increases; an explanation of this 
statement was given later by S. Timoshenko (2). The representation of 
th deflexion-curve in the form of trigonometric series is especially useful 
in this case and leads, as is shown below, to the very simple formula (6) 
giving the general solution of the problem for any number of spans. 


Details of numerical solution 
We start from the known differential equation for the deflexion-curve 
of a beam submitted to the simultaneous action of a lateral load and of an 
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BUCKLING OF CONTINUOUS BEAMS ON ELASTIC SUPPORTS 259 
axial compressive force P (Fig. 5): 
Elu"+M = —Pu. (1) 


Here EJ denotes the flexural rigidity of the beam and .@ the bending 
moment produced by the lateral force Q alone. It can be represented in 


the form of the trigonometric series 


Ma @ p. asin, sin“ 











or, if we replace Q by —u,/2 (u; being the deflexion at x = a,): 


2i —. 1 . kna; . kax 
M ~~ S <a, 
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We assume that the deflexion-curve can also be represented by a 


trigonometric series : 
. . nn 
. . ° ) 
‘ S Cn sin 1 ? (2) 
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n 1 
: 3 5 . _nmewl. kra; krx 
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The bending moment produced by the simultaneous action of the lateral 
forces R,, R,,..., R,,-, will be 


-_—" a 
ae 2s 4 “ . na, <1. kra; koa 
M > A : p aR.” c, sin : > - sin ——_sin ——, 
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Introducing the expressions (2) and (3) into (1) 
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co m—1 
P és . ka; . nna; P 
we get n'c,, +2Q > Ch > sin——sin—— = yn*ey, (4) 
k=1 i=1 


where the following notations are introduced, 


PI? 
QO = —-, ? = a 
mEIU Cd Oe | 
If the intermediate supports are equally spaced, a; = il/m, and it 
follows from (4) that 


8 
nic, +Q > ce, K(n,k) = yn’c,, (5) 
k=1 
m—1 e ° 
: . nin . kin 
where Kin, k) = 2 > sin ——__ cia —. 
— m m 
v 


This coefficient vanishes when k or n are integral multiples of m. If they 
are not, then 


m—1 


K(n,k) = > = 


u(k—n ar a i(k nmr, 
| m m | 


i=] 
Using the known expression for trigonometric polynomials 
mal ; 1 1 sin(m—4)0 
> cosid = —~4+-=—____-_, 
‘<1 2° 2 sin3é 
it is easy to prove that K(n,k) = m, if k = n+2jm (7 = 1, 2,...) and that 
K(n,k) = —m,ifk = —n+2jm. For all other values of & the coefficient K 
vanishes. 
Thus for values of x which are integral multiples of m we obtain from (5) 
a group of independent equations 


Bn? = y, 


22m? = y, 


corresponding to the buckling of the beam into m, 2m,... half-waves. The 
smallest value of the critical force P,, is given by the first of them 


m 
Yn c= m*, 
which corresponds to the buckling into m half-waves. 

For values of x which are not integral multiples of m we obtain homoge- 
neous linear equations in the coefficients c,,. These equations can be divided 
into m—1 independent systems. Setting their determinants equal to zero 
we obtain m—1 values of the critical force, corresponding to buckling into 





We 
force 


[(m— 


We 
colun 
of th 


mo 


or 


wher 


they 


that 
it K 
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1, 2,..., m—1 half-waves. For example, the equation corresponding to the 
buckling into m—r half-waves would be 


n—r)'§—(m—r)*y+mQ mi mQ  % 
mQ (m+r)*—(m+r)}y+mQ —mQ . e 
mQ mQ (3m—r)*—(3m—r)*y+mQ. .|= 0. 





We impose the condition that this force must not be smaller than the 
force required for buckling into m half-waves and replace y by m? 


m—r)?—m?|(m—r)? + mQ mid mQ ‘ 
mQ [(m+r)?—m?}(m+r)?+mQ — mi . « 
mQ mQ [(8m—r)?—m?](8m—r)?+mQ. .|= 0. 





We retain 2x columns and 2n rows in the determinant and add every 
column to the foregoing one. Then, developing the determinant in terms 
of the constituents of the last column and increasing n to 00, we obtain 


m2 5 } . a —>( — —l, 
LH, \{(sm—r)?—m?](sm—r)? © [(sm-+-r)?—m?](sm-+-r)? 
s=1,3,5 , 
(5°) 
4m?—r?)(m?—r?)? , 
or Q(] R) — _ Mi as ) P (6) 
2m(5m*-+-r*) 
where 
R (4m?2—r?)(m?2—r?)* — | ] 


2(5m?-+-r?) Lz || (sm —r)®@—m?](sm—r)? 
s=3 
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+ : re cy 7 oO ) ! wT) a“ 
[ (sm-+-r)?—m?]|(sm-+-r)? 
The greatest value of R corresponds to r 1, and 
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262 BUCKLING OF CONTINUOUS BEAMS ON ELASTIC SUPPORTS 
so that it can be omitted from (6). Thus the greatest value of Q is 
(4m?2— 1)(m?—1)? 





0 om(om*1) (6) 
73 
and a — 6" 
Instead of UX, I. G. Boobnov used 
iI m4(5m?-+- 1) 


9 

ss en peed ee 
0 = mime t= 38 (mt —1)(ame—1)" 
Conclusion 

The comparative table appended below gives the results of equation (7) 
and the values of ® as calculated by I. G. Boobnov. Increasing m to 
we obtain from equation (7) the value 0-253; but if we also retain the 
second term of the series (5’), the value 0-250 is obtained, which agrees 
with the value indicated by I. G. Boobnov. 

















m= | 2 3 | 4 | 5 | 6 7 9 II 
I. G. Boobnov 0500 | 0°333 | 




















. . | 
Equation (7) | 0°504 | 0°337 | 0296 | 0280 | 0-271 | 0°265 | 0261 | 0°259 
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ANALYSIS OF CONTINUOUS STRUCTURES BY THE 
STIFFNESS FACTORS METHOD 
By LEROY A. BEAUFOY and A. F. 8S. DIWAN 
(College of Estate Management, 11 Great George Street, London, S.W.1) 
[Received 16 October 1948] 
SUMMARY 


The paper presents an ‘exact’ method of analysis of highly redundant two- 


limensional continuous structures. No simultaneous equations are used, and all 


calculations may be done by slide rule. The method is a mathematical analogue 


of the experimental one using unloaded models. The general case is considered, 
but for ordinary use considerable simplification is possible. The procedure is 
illustrated by an example of a three-span continuous-arch system involving nine 


redundancies. 


Introduction 

THE stiffness factors method was developed using the concept of a hypo- 
thetical unloaded model imagined to be manipulated in the same way as 
an actual model (1); influence lines may be found, or moments and forces 
for a given loading computed directly. The distinctive characteristic of 
the experimental method is retained, that is, a structural analysis is first 
made of the unloaded structure. This analysis, which yields information 
about the deformation of the structure as a whole due to any given 
imposed distortion, is viewed as an expression of the elastic properties 
of the structure, which are independent of the loading carried by it. Once 
these are known the effect of any loading which does not strain the 
structure beyond the elastic range is easily investigated. Joint displace- 
ments in rotation and translation are found. For a loaded structure, 
these are obtained as direct displacements due to unbalanced fixed-end 
reactions on a joint plus displacements induced at the joint by balancing 
other joints. The total joint displacement being known, end moments 
and forces acting on members can be computed without using simultaneous 
equations or successive approximations. 


Definitions 

Stiffness factors 

The stiffness factors for a member end or joint are the moment and the 
horizontal and vertical forces required at that end, or joint, to produce 
unit displacement there. These factors are called horizontal displacement, 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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vertical displacement, or rotation stiffness factors according to the direction 
of the displacement. 

The stiffness factors for the end A of a member AB in a structure are 
called restrained or actual respectively when the far end B is assumed 
fixed or free to move (together with all other joints in the structure on 
the same side of A as B) to the extent which the stiffness of the members 
permits. 

Stiffness factors for joint B are called modified or actual respectively 
when they give the stiffness of joint B when joint A is assumed fixed or 
free to move in conformity with the stiffness characteristics of the struc- 
ture as a whole. 


Displacement coefficients 

The displacement coefficients for a member end or joint are the displace- 
ments due to unit moment or unit horizontal or vertical force applied at 
the end or joint. The displacement coefficients for joint B are called 
modified or actual respectively when the neighbouring joint A is assumed 
fixed or free to move in conformity with the stiffness characteristics of the 
structure as a whole. 

Induced movements 

A unit rotational, horizontal, or vertical displacement applied at joint 4 
will, in general, induce a movement of joint B, which can be resolved into 
rotational, horizontal, and vertical components; these are called induced 
movements. In turn, they induce further movements at the next joint C, 
and so on. 


Notation and sign convention 
Symbols used without subscripts 
¢@ rotational displacements 
A horizontal displacements 
A vertical displacements 
m, M moment 
h, H horizontal force 
v, V vertical force 
2, y horizontal and vertical coordinate distances, respectively 
&, 9 horizontal and vertical coordinates respectively of the elastic centre 
El flexural rigidity of a structural member or unit 
S elastic area of a structural member or unit 
ds elementary length of a member 
AB member AB; end A of member AB 
a determinant 
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ANALYSIS OF CONTINUOUS STRUCTURES 
Subscripts 
¢ denotes that the quantity is in respect of applied rotational dis- 
placements 
A denotes that the quantity is in respect of applied horizontal dis- 


placements 


> 


denotes that the quantity is in respect of applied vertical displace- 
ments 


” 


denotes that the quantity is in respect of applied moment 


= 


denotes that the quantity is in respect of applied horizontal force 
» denotes that the quantity is in respect of applied vertical force 


Supe rscripts 
be denotes that the quantity is in respect of the member BC 
3 denotes that the quantity is in respect of the joint B 


Symbols used with subscripts d, A, A 
m,h,v Moment, horizontal force, and vertical force respectively, 
acting simultaneously on one end of a member to produce 
unit displacement there when the other end of the member 
is restrained in the same manner as in the structure (actual 
stiffness factors for member ends) 

m',h',v’ As for m, h, v but when the other end of the member is 
considered fixed (restrained stiffness factors for member 
ends) 

M,H,V Moment, horizontal force, and vertical force respectively, 
acting simultaneously at a joint to produce unit displace- 
ment there when all the other joints are restrained in the 
same manner as in the structure (actual stiffness factors for 
joints) 


M’',H’', V’ As for M, H, V but when the far end of one of the members 


connected to the joint is considered fixed (modified stiffness 
factors for joints) 

¢, A, A Rotational, horizontal, and vertical displacement respectively, 
produced at one joint of a structure by unit displacement 
applied elsewhere in the structure (induced movements) 


Symbols used with subscripts m. h.v 


¢, A, A Rotational, horizontal, and vertical displacement respectively, 
produced at one joint of a structure by unit moment or force 

applied at the joint (actual displacement coefficients) 
¢', \’,’ (Modified displacement coefficients corresponding to the 


above) 
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Others 
e Carry-over factor for moments or forces 


=~ 


, Centroidal moment of inertia of a member or structural unit 
about the x-axis 

I, Centroidal moment of inertia of a member or structural unit 
about the y-axis 

I, Centroidal product of inertia about axes a, y 

I = L—(i3 [,) 

I, 


xy 
= L,—(Iz I.) 


xy 
Sign convention 
Moments and rotational displacements are taken as positive when 

clockwise; horizontal and vertical forces, displacements, and coordinate 

distances as positive when measured to the right and upwards respectively, 


Displacements of an elastic member 


Let Fig. 1 represent a member AB, of flexural rigidity EJ, fixed at B. 





Fie. 1. 


Imagine a rigid arm to extend from end A to the elastic centre O; at 0 
let there be a vertical force v, a horizontal force h, and a moment m; 
take O as the origin of axes Ox and Oy. 

When m is unity and acts alone we have 


$= [ds/EI=S (1a) 

A=—|yds/EI =0 (16) 

A: [ xds/EI = 0. (1¢) 
When A is unity and acts det we have 

¢= — [ yds/EI = 0 (2a) 

A = [ y*ds/EI = I, (2b) 


A= — [ ay ds/EI = —I,,. (2c) 
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When v is unity and acts alone we have 


¢= | xds/EI = 0 (3a) 
A { ay ds/EI = —L,, (3b) 
d= | a ds/EI = I, (3c) 


then m, h, and v are given any values and act si aneously the tote 
WI h, and v are giv ; | nd act simultaneously the total 


displacements of the elastic centre are given by 


db m S (4 a) 
A = hI,—vl,, (4b) 
A vl, h iw (4c) 


Transfer the origin of the x- and y-axes to end A and let the coordinates 
of the elastic centre O be % and 7. Then the displacements of end A are 


given by d ms (5a) 
A hl,.—vI,,—Y9 (5b) 
A vl, —hIL,, + $e. (5c) 


The quantities S, J, [,, and I, above, together with the coordinates Z, 7 
of the elastic centre with respect to end A, are called the elastic constants. 


Restrained stiffness factors 
Suppose (Fig. 2) that a moment m and forces h and v are required at O 





Fic. 2. 


to produce a unit rotation of end A without translation, that point O 
moves to O,, and that the horizontal translation of O is 9, and the vertical 


translation is (—#). From equations (4) we have 
d ms l (6a) 
A hI vly q (65) 


A vl, h = —Z, (6c) 
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from which we deduce that 





ae WE (7a 
y tL y _ 
h = — lr] (7 h 
ay Fo ‘, 
UL x y 
Y == Se, (7c 
fate By 
The corresponding values at A are given by 
"2 72 One 
m= +545 oe (80 
Ss I; ly l;, [, 
j $i. 
hy = in (8b) 
fa Sale 
, GLy A 
uw = SS TF (S¢} 
9 « iy 


These values are the restrained rotation stiffness factors for end A of 
member AB. 

Horizontal displacement 

Suppose that a unit horizontal displacement is imposed on end 4 
without permitting rotation or vertical translation there; the moment 
m, and forces hi, and v, required at A, viz. the restrained horizontal 
displacement stiffness factors there, may be shown to be 


, y 7. ( 
Mx (9a) 
: % ie 
Ron iff, (9b) 
v, = L,,/I1,1,. (9c) 


Vertical displacement 


Similarly the restrained vertical displacement stiffness factors are given 


by “ -. 
‘ ' iL z 
m = a_—— (10a) 
. iy i. I, 
hy = Ly I; I, (106) 
v, = 1/J,. (10¢) 


As the far end B is fixed, the restrained stiffness factors given above are 
also the actual stiffness factors for end A. 

Horizontal prismatic member 

For a straight horizontal prismatic member A B of length L and constant 
flexural rigidity HJ, the elastic constants and stiffness factors for three 
particular cases are given in Tables 1 and 2 respectively. 
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TABLE 1. Elastic constants for a horizontal prismatic member 





Condition of right-hand end 


Elastic Freely 
constant Fixed Hinged supported 
S L/EI we o 
i L/2 L/2\1 0|;L 0 
I 0 0 CO 
I L3/12EI L3/3EI L3/3EI 
I 0 0 0 


TABLE 2. Stiffness factor values for left-hand end of a horizontal 
prismatic member 


Condition of right-hand end 


jactors Jo? 
left-hand Freely 
end Fixed Hinged supported 
m’ $I /L 3E1I/L 3EI/L 
h’, = m’, 0 0 0 
; m, 6E£I/L* 3EI/L? —3EI/L? 
/ : x 4 0 
h’ — v’ 0 0 0 


12EI/L 3E£1/L* 3E1/L* 
Joint stiffness factors and balancing equations 

The actual stiffness factors for a joint are the sums of the relevant 
stiffness factors for the ends of all members connected with the joint, 
since any unit displacement imposed on the joint also imposes unit 
displacements on all such ends. Thus, for example, the actual rotation 
stiffness factors are 

M, = > m,; Hy = > hg: Vg = > 0%. 

To determine the induced movements at the joint B and the actual stiff- 
ness factors for the end A of member A B (Fig. 4), modified stiffness factors 
for the joint B, which give the stiffness of that joint when the joint A is 
wssumed fixed, are used. They are equal to the sums of the actual factors 
ior the ends of all members connected to the joint B, excluding BA, plus 
the restrained factor for BA. Thus 


M;, my mil } mise (lla) 
H', = WE +hg +h. (115) 


The modified stiffness factors for a joint therefore depend on whether it 
is approached from the left or the right. 
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Since the actual stiffness factors for a joint represent the moment and 
forces required to produce a unit displacement there, it follows that if the 
joint has simultaneous movements ¢, A, and A, the required values M. # 
and V are given by the following balancing equations 


M = $M3+AM,+AM, (12a) 
H = $¢H,+AH,+AM (126) 
V = ¥+A+Ah. (12¢ 


Displacement coefficients 
When M is unity and acts alone we have 


$m My+An My +A, My = 1 (134) 
4. ig +-A,, H+, Bh = © (130 
é.%+4,,%.+4, KR = 6. (130) 
If k=|M, MM 
H A, Ay 
‘Rn ws 


= My(H, \—H})—Hg( Hg W—Ay Vy)+-V( Hy Ay— HV), 
in which use is made of the reciprocal relations M, = Hy, M\ = V;,, V, = Ay, 
the solution of equations (13), representing the displacement coefficients 
for moment, can be written 


dm = (Ayn -H})/k (14a) 
4, F — (Hs V,— Hy Vz) k (145) 
An = (Hy Hy—HyV,)/k. (140) 


These values, together with those for the displacement coefficients for 
vertical and horizontal forces, which can be found in a similar manner, 
are collected together in Table 3 (column 2). Expressions for special cases 


TABLE 3. Displacement coefficients in terms of stiffness factors 


Nature of joint displacement 


2otation and Rotation Rotation 


Displacement translation in any and vertical and horizontal | Rotation 

coefficient direction translation only| translation only| only 
(1) 2) (3) (4) (5) 
dm 4 /(AM;,—BH,-+-CV4) V,/D H,/E 1/Ms 

An = C/(AM;,—BH,-+-CV) —V/D 0 0 

A, = ody — B/(AM,—BH,+-CV;) 0 H,/E 0 
Ap D/(AM,—BH3+CV;) 0 M,/E 0 
Ay E/(AM,—BH,-+- CV) M3/D 0 0 

A Ap, F/(AM,—BH,+CYV,) 0 0 0 


,—-M; B=HVy—-Hh; C=HH,—Hh; D=M%y%-" 
E=M,H,—Hi; F = Hy V%—MgH, 


—_——— 
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of restricted movements of joints are easily obtained from them (columns 
3-5), If actual stiffness factors are used in this table, the resulting displace- 
ment coefficients give the actual displacements of the joints; if modified 
stiffness factors are used, the resulting displacement coefficients give the 
displacements of the joint for the condition to which the modified stiffness 
factors apply. 

When MV, H, and V are given any values and act simultaneously, the 


halancing movements at the joint are 


¢ = M¢,,+H¢,+V4, (15a) 
A = MA,,+HA,+VA, (155) 
A= MA,,+HA,+D),. (15) 


Actual stiffness factors and induced movements 

When the end B of the member AB is not fixed, but is connected to 
one or more members for which the actual stiffness factors are known, 
the actual stiffness factors for the end A are obtained as follows: First 
fix the end B and impose the desired unit displacement at the end A. 
Next, fix the end A in this position and release joint B so that it takes 
up a position of equilibrium determined by the modified stiffness charac- 
teristics of joint B with respect to the end A. Then the forces and moment 
acting at the end A are the actual stiffness factors for AB. 


Rotation 

When the end B is fixed (Fig. 3) a moment mz and forces hy and 4 
are required at the end A to cause unit rotation there. These set up an 
unbalanced condition at B, due to a moment m and forces h and v such 


that 
. ‘ ab ’ ab / ab» > 
m moa —hg™y, +0372, (16a) 
h hie (165) 
v yr, (16c) 
Further, m c%m',, in which ¢ is the carry-over factor for moment 
trom end A to end B due to rotation at A. (For prismatic beams, c 4.) 


Lock the joint A in this position and release joint B; the balancing move- 
ments there are the induced movements at B due to unit rotation imposed 
it A, These induced movements are found by substituting the values 


obtained from equations (16) in equations (15), giving 


D4 ( omg din | hig’dy, | v6" >, (17a) 
Ag ee mie’, T hg?A,, | vg™A, (175) 


—bmi,%X, | her, | vg" X,. (17 Cc) 

















272 L. A. BEAUFOY AND A. F. S. DIWAN 





These movements at B require the action on member BA of a moment m whi 


and forces h and v given by 
m = ogmg*+Ag miyo4-+-Ag myo4 
h = dg hgrt+Ag hinot+rAg hy 
v= bg vg +Ag vyt+Ag v4. 








~< xy ae ee ee ge 
| / 
/ 
| 4 
| ae 
Pat 
ae 
~ 
~ hi} 
iy oe — 
at 
,ad 
Vv 
Y Fic. 3. 


These actions at B produce reactions at A given by 


m = cd. m40¢— Ag m2 —Ag my” 


h _ ds he’ -Ag hiw— rz hye 
v —dbzy vyea a Ag v,™ = rs v7, 


Hence, the total forces at A which, since B has now been released and 4 
has been rotated through unit angle, are the actual stiffness factors for 


the end A of member AB, are given by 


/ ’ ba , 
/ 


, 
; m m m 
M4, - mI +chag, —¢ : —A,—a—A =) 


mi, mis om f 


QD Cea) 
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hg ng(1— gt Maze) 


he Ph PR 
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Us vs(1— dg #_—A,-A_A,-). 
“\ U} "U4 US 


Horizontal displace ment 


Similarly, the induced movements at B due to a unit imposed horizontal 


displacement at A are 
dy md’, | hid, | vd, 
bad! PabA’ 1 a! ab? 
A, mM ~< = : h* An rv," A 


, ba , , ah , ,/ uh, , 
Mr, thr, + vA,» 
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while the actual stiffness factors for end A are 





af m m 
ma m'{ 1+ c'"g,, = —Aa— rx m) 
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: jibe h‘ 
hs his(1— ry —A,—-Ay 
TIN YA 
; v ' ba v5 
vy = val 1—d, S$ — Ady), 
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Vertical displacement 
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The induced movements at B due to a unit imposed vertical displace- 


ment at A are tha J? 

nent dy myadg’ hyd,,4 -vydi, 
Ay = myPaA’, +hyPAy + ope, 
A) my"*r,4 hy, ; v,X,,, 


while the actual stiffness factors for end A are 


f ‘ba J 
, m m 
my = m){1 +0", —£—— Ay— > — ») 

my my 
hi,be h', 

hy | l o,— nas Az - Ay 
hy “hs 
A A 


Axial deformations 
When axial deformation of 
member is connected to a joint, A 


equations (17) and (21), are 


dd om gaig’ hig’ dy dy 
Ag 7m A. his’) 
by = My*h, hn dy, 

Ay mA’, . hi, Ay, 


members is disregarded, 
, the induced movements, 


(23a) 
(23 b) 
(23c) 


(24a) 


(245) 


(24c) 


and a vertical 
from 


in which the displacement coefficients in terms of the modified stiffness 


factors are as given in Table 3, column 4. 


The actual stiffness factors 


m4 


, from equations (20) and (22), are 
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Symmetrical members 

For symmetrical arches and bents, when the restrained stiffness factors 
for the two terminals are equal (except in the cases of m‘, and v4, when the 
signs will be reversed), if the span A B is written as LZ and the height of the 
centroid O above the (horizontal) springings AB as 7, then 


—(m4+-v4L) 








cab — : ale 
ms 
hi l 
his 2 i 
and equations (26) reduce to the following: 
, ma 5 
m m +-cd3— OT 
¢ (1 $ 2) (2/4) 
A = 
hg = Wis 1—45— =!) 27h 
Y 
m’ 
ma ms (1 + ch, “¢—A,] (27¢ 
my 
hy = hy(l1—¢59—A,). 27d 


Total end reactions on a member 

For a member AB subject to imposed distortion let the displacements 
at A be ¢4, A“, A4, and those at B be 48, A®?, X*. Then the moment 
and forces required at A are 


M4 = 64mg” +8 c%mgot + (A4—A)mi, + (AA—AP)my” ~— (28a 


H4 = $4hia— pb hgot4 (A4—AB)hAe+ (A4—)B)Hje (285 
VA d4y,% _fP vot | (A4- AP )v', f° ab_t (A 4 __ By \a0. (28¢ 
while those at B are 

MP2 pPmyet | d tommy” | (A2 A {)m/,°* it (A¥ A+ )m)e2 (29a 
He hBh',ba— pAp',2v 1 (AB A- ty pia 1 (\B__}A)fjba HA (29b 

d db A 
VB hBy',ba_ fy jab AB A\A)py',ba_1. (AB )A)yiba VA, (29¢ 

¢ Up r 


[f, in addition to imposed distortion, the member AB carries external 
loading, the fixed-end reactions due to this loading must be added to the 
above values. 
Description of the method 

Direct analysis consists of three distinct steps as follows: 

Step 1 

Divide the structure into units which are first considered as fixed- 
ended. Find their elastic constants and restrained stiffness factors. Then 
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determine the actual stiffness factors and induced movements of the joints. 
The structure is considered unloaded, but subject to the restraint con- 
ditions obtaining at its supports, so that its deformation depends only on 
the elastic properties of the structure itself. This step constitutes the 
structural analysis. 

Step 2 

Assuming all joints locked, determine the unbalanced moments and 

forces at the joints due to the given loading. Then release all joints and 
find the balancing movements at each one in turn due to its own previously 
unbalanced moments and forces (the other joints being unrestrained), 
together with the resulting induced movements at other joints. The total 
displacements at any joint are then the sum of balancing movements at the 
joint plus movements induced there by balancing movements elsewhere. 

Step 3 

Compute moments and forces throughout the structure as follows: 

Moment = (m4, total rotation at joint)+(m) x total horizontal dis- 
placement at joint)+-moments carried over from the far end of the 
member due to rotation and horizontal translation there. 

Horizontal force = (hy x total rotation at joint)+(h4 x total horizontal 
displacement at joint)+horizontal forces carried over from the far 
end of the member due to rotation and horizontal translation there. 

These expressions follow from equations (28) and (29) and refer to the 
case where there is no vertical translation of the joint. 

Steps 2 and 3 together constitute the process of stress analysis. 


Example 

The method will be applied to the continuous-arch system shown in 
Fig. 4 which has been solved elsewhere (2) by the method of successive 
approximations. The system is composed of three similar elliptical arches 
of 60-ft. span and 15-ft. central rise, carried on piers 20 ft. high. The 
arches and the piers are of the same constant cross-section. 

\ssume the curved length of the elliptical arch to be divided into 
eleven equal lengths each of 6-8 ft., and the coordinates of the mid-points 
of these lengths to be as follows: 





3070 | 36°55 | 43°10 | 49°35 | 55°10 | 59°30 ft. 
: S| & | 59°; 


15°0 | 14°6 13°5 | II's 8-2 | 3:2 ft. 





the distributions of moments and of horizontal forces throughout the 
structure are to be found for a uniformly-distributed loading of 1,000 lb. 
per ft. run covering one of the end spans, neglecting axial deformation. 
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The complete solution is given in Tables 4, 5, and 6, to which the fo]. 
lowing explanation refers. 

Ste Pp 1 

Calculate the elastic constants for a single arch. Since ds/EI for the 
eleven equal lengths referred to above is 6-8, if unit value is assumed for 














» — aces 3 BE GOR: cach = OK, > 


Fic. 4. 


EI, then S = 74:8, the total curved length of the arch. The coordinates 
of the elastic centre are = +30 ft. by inspection, and 

j= >dy/1) 10-63 ft. 
Transfer the a-axis to pass through the elastic centre and obtain coordi- 
nates corresponding to those given above. Then J, = 6-8 ¥ y? = 1,300; 


I 6-8 ¥ 22 = 28.300: J 
—— 


y -y = 9, since the arch is a symmetrical section. 


The elastic constants are now known for the arch and are entered in 


Table 4, item (1), columns CD and DC. For any pier such as DH, the 


elastic area = L/EI = 20; & by inspection = 0, and 7 10; 
IT, = [°/12EI = 667; 

[, and I,,, are both equal to 0; these values appear in item (1), column DH 
of the table. Using equations (8) and (9), obtain values for the restrained 
stiffness factors for the separate arches and piers and record them in 
item (2). The flexural rigidity of the arch and pier sections was taken 
as 1 kip-ft.2, which keeps the units consistent, but does not otherwise 
affect the solution. 

Starting at the right-hand end of the structure, obtain by summation 


the modified stiffness factors for joint D as follows: 
Mi, = 0:20+-0-1324 = 0-3324 
Hy = —90-015-+-0-0082 = —0-0068 
Vg = 0-0015+-0-00077 = 0-00227, 


and enter these values in the table in item (3), column DH. Next calculate 
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the modified displacement coefficients for joint D from Table 3, column 4 
substituting the modified stiffness factors, and enter them in Table 4 
item (4), column DH. 

Now calculate, from equations (25 a and 6), the rotation 44 and transla- 
tion A, induced at D by imposing unit positive rotation at C. Thus, 
dy = 0°30; Ag = 4-52. Enter these values in the table in item (5), column 
DH, (The value of the carry-over factor c for the arch, to be used in 
these equations, is obtained as follows: Consider the arch CD with a 
fixed end at D. A moment m4, and forces hg, Us are required at C to produce 
unit rotation there; a moment m” is set up at D. The value of mg = 0-132 
is known (Table 4, item (2)), and 7% = —0-001063 can be found from 
equation (8c). Then moments about D give 

mP = 0-1324—0-001063 x 60 = 0-0684. 
Hence ¢ = —m?/m = —0-0684/0-01324 = —0-517). 

Similarly, calculate from equations (25 c and d) the movements ¢, and 
A, induced at D by imposing unit positive horizontal translation at (. 
Thus ¢, = 0-034; A, = 0-44. Enter these values in the table under 
item (5), column DH. Now derive values of the actual stiffness factors 
for CD from equations (27 a, b, d) as follows: 

ms = 0°1324[1—0-517 x 0-30—4-52(0-0082/0-1324)] = 0-075 

hg = 9-0082[1—0-30—4-52(0-00077/0-0082)] = 0-00227 

hy = 0-00077[1—0-44—0-034(0-0082/0-00077)| = 0-000154, 
Enter these values in the table in item (6), column CD. 

Note that the actual stiffness factors for DH, CG, BF, AE are the 
same as the restrained stiffness factors since in each case the far end of 
the member is fixed. 

Next calculate the modified stiffness factors for joint C. They are: 


My = mF +-mF +m = 0-075+0-20+0-1324 = 0-4074 
Hy = hy th? +hg = 0-00227—0-015+-0-0082 = —0-00453 
Hi, = h@+-hO+hie” = 0-000154+-0-0015+-0-00077 = 0-002425. 


Enter these values in the table in item (3), column CG. 

The cycle of operations performed above may be described as follows: 
Starting with the modified stiffness factors for joint D, find the modified 
displacement coefficients and then the induced movements at the joints 
and hence the actual stiffness factors for CD, which enable the modified 
stiffness factors for joint C to be found. This cycle is repeated, giving in 
turn the modified displacement coefficients and the induced movements 
at joint C, the actual stiffness factors for BC, and the modified stiffness 
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factors for joint B. A final cycle gives the modified displacement coefti- 
ients and induced movements at joint B and the actual stiffness factors 
for AB. This enables items (3), (4), (5), and (6) in the table to be completed 
ig shown. 

The actual stiffness factors for the joints, item (7), follow by summation 
of values in item (6), and all required stiffness factors are known. Item (8) 
gives the actual displacement coefficients, which are obtained from Table 3 
by substituting in it the actual stiffness factors. 

Next prepare the unit displacement table, item (9), which shows the 
lisplacements ¢@ and A at all joints due to unit imposed rotation or 
translation at any joint. This part of the analysis is now complete, and 
t may be used to determine the effect of any particular loading. 

Step 2 

Consider the uniformly-distributed loading of 1,000 lb. per horizontal ft. 
m the end span AB. The fixed-end reactions are: 

Fixed-end moments (ft.-kips), 126 at AB; —126 at BA 
Fixed-end thrusts (kips), 36 at AB; —36 at BA. 


Elsewhere they are zero (see Table 5, item (1)). Hence, the unbalanced 


components at the joints are: 
For joint A: Moment —126; thrust —36 
For joint B: Moment 126; thrust 36. 


These values constitute item (2) in the table. 
Find the balancing movements at the joints from equations (15 a and 6). 

Thus, for joint A: 

d 126 x 4:-44—36 x 26-7 = —1525 

A 126 x 26-7 —36 x 705 = —28850. 
Enter these values in Table 5, item (3). Next use the unit displacement 
table to obtain by proportion the movements induced at other joints by 
these balancing movements at joint A and enter the results in item (3) in 
the table. Apply the same process to the unbalanced components at 
joint B, and enter the values as item (4). Sum the values in items (3) 
and (4) to obtain total joint displacements, item (5). 


Step 3 

The total joint displacements known, computation of moments and 
thrusts may be made, using equations (28). This is done in Table 6, which 
will not require further explanation. It will be seen that agreement with 
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Maugh’s results is very close indeed, although calculations were made 
only with a 10-in. slide rule. 


Advantages of the method 

Principal among these are: 

1. The total elimination of simultaneous equations, coupled with the 
fact that the method is an exact one; as a result, computations may be 
performed throughout by slide rule. 

2. A single structural analysis, independent of the loading, is first made, 
From this, the effects of any loading are readily found. 

3. Translation, as well as rotation, of the joints is considered, so that 
the effect of side sway is automatically taken into account. 


4. Frames containing curved members of variable cross-section are 
readily analysed. 
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THE INDUCTION OF ELECTRIC CURRENTS IN 
NON-UNIFORM THIN SHEETS AND SHELLS 
By A. T. PRICE.(Imperial College, London) 

[Received 12 October 1948] 


SUMMARY 

General equations are obtained for the induction of electric currents in any thin- 
sheet distribution of conducting material by periodic or aperiodic fields. Methods 
of solving the equations for non-uniform plane sheets and spherical shells are 
considered. Numerical calculations for a special case of the plane sheet are made 
to illustrate some important features. The case of the spherical shell is considered 
in some detail, with a view to geomagnetic applications. Methods of approximating 
to the solution when the conductivity is given as an empiric function are described. 


1, Introduction 

Tae theory of the induction of electric currents in non-uniform thin-sheet 
distributions of conducting material is of importance in several geo- 
magnetic applications. When, for example, attempts are made to infer 
the distribution of electrical conductivity within the earth from analyses 
of the observed daily variations or storm-time variations of the earth’s 
magnetic field, on the assumption that the part of the varying field of 
internal origin is due to currents induced in the earth by the part of 
external origin, it is found that the oceans, which are relatively good 
conductors, must have some appreciable influence on the internal part 
Chapman and Price, 1, 1930; Lahiri and Price, 2, 1939). To determine 
the exact nature of this influence it is desirable to calculate the currents 
induced in a thin spherical shell having a suitable non-uniform distribution 
of conductivity. The general theory of thin spherical shells is dealt with 
in § 9 of this paper, and the analytic solution for a special case is obtained 
in §§ 10 and 11. Methods for approximating to the solution for other 
distributions of conductivity are described in §§ 12 and 13. 

Similar calculations relating to non-uniform plane distributions would 
be useful in investigations relating to earth currents, including those 
flowing in sea channels (Barber, 3, 1948), when account is taken of the 
local variations in resistivity of the earth’s crust. The theory of plane 
sheets is dealt with in §§ 4 to 6, and numerical solutions which illustrate 
some important features are obtained in §§ 7 and 8. 

Another important application of the theory is the investigation of the 
effects of the non-uniform distribution of conductivity in the ionosphere 
on various geomagnetic phenomena, including sudden commencements, 
micropulsations, and short-period fluctuations of the earth’s field during 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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magnetic storms. These have recently been discussed in a joint paper by 
A. A. Ashour and the writer (4, 1948), and the numerical calculations 
therein described are based on the formulae obtained in §§ 9 and 10 of the 
present paper. 


2. Statement of the general problem 

We consider any non-uniform thin-sheet distribution of conducting 
material, the medium on both sides being non-conducting and the magneti 
permeability everywhere unity. A varying magnetic field in the surround- 
ing medium, having its origin in sources external to the sheet, e.g. produced 
by electric currents in a region not in contact with the sheet, induces 
electric currents in the sheet. It is required to find the magnitude and 
distribution of these currents, and also their magnetic field. 

We assume that the variations of the field are sufficiently slow to enable 
us to ignore the displacement current and treat the velocity of propagation 
of electromagnetic effects in the dielectric as infinite. The permissible 
rapidity of variation of the field will then be limited by the average con- 
ductivity of the conductor and by its linear dimensions. 


3. General equations for any thin-sheet distribution of conducting 
material 
Let H, E, and i denote the magnetic intensity, electric intensity, and 
the conduction current density, respectively, all measured in electro- 
magnetic units. Neglecting the displacement current, the conduction 
current is non-divergent and the field equations therefore reduce to 


curlH = 4zi, (1 
curlE = —p éH/ét, (2 
i= «EK, (3 


where « is the specific conductivity, and is a function of the space variables. 
The magnetic permeability u is assumed to be unity everywhere, so that 
div H is zero. In the dielectric curl H is zero, hence 


H = —gradQ, (4 


where the potential Q satisfies Laplace’s equation. 

On passing from the dielectric into the conducting material, the normal 
component H,, of H is continuous, since H is everywhere non-divergent 
It follows that, if the thickness d is infinitesimal, the normal componett 
H,, will have the same value at corresponding points on opposite sides of 
the sheet. 

By applying the circuital relation (2) to a small rectangular circut, 
having two sides parallel to the surface of the conductor, one being in the 
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conductor and one in the dielectric, it is seen that the tangential compo- 
nents of E are continuous at the interface. Hence, for an infinitesimally 
thin conducting sheet the component E, of E, parallel to the sheet, has 
the same magnitude and direction at corresponding points on opposite 
sides of the sheet. 

Integrating equation (3) over the thickness d of the sheet at any point, 


we obtain d d 
| idZ | x«E df, 
0 0 
q 
where the integral on the right may be replaced by E, | « df to the first 
0 


d 
order of small quantities, when d is small. The integrals | id{ and 
0 
q 
| «df represent, respectively, the total sheet-current intensity per unit 
) 


length of surface at the point, and the total conductivity per unit length; 
denoting these by i, and 1/p respectively, we have 

E. pi,. (9) 
Applying equation (2) to a small circuit in the surface of the sheet, we 


ob re i ~ 9 lage 
— curl E, n 0H, /et = n &?Q/eten, 


where n is a unit vector normal to the surface from the negative to the 
positive side. Substituting for E, from (5), this becomes 
p curli,+- grad p / i, n 0H, /et = n €?Q/eten. (6) 

If we now apply the circuital relation (1) to a small rectangular circuit, 
having two sides of length 5s parallel to the conducting sheet, with one 
on each side of it, we obtain 

(H,).,. 88—(H,)_ 8s = 471, ds sin 8, 

where the + and — suffixes indicate values at the positive and negative 
faces of the sheet, and @ is the angle between ds and i,. Since this relation 
is true for any direction of 5s tangential to the sheet, we can write 
os = 
i 0 (H,—H_) = —nA grad, (7) 
— , l 
where y —(Q.,—Q_). (8) 


Equation (7) shows that ¥ is the stream-line function for the currents 
flowing in the sheet. 
Substituting the value of i, from (7) in equation (6), we get 
pcurl{n \ (H.,—H_)}+ grad pA {n, (H,—H_)} = —4an ¢H,,/ét, 
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which is easily reduced to the form 

pdiv(H,—H_)+ grad p.(H,—H_) = —4ré H,, /et. (9 
This equation gives the relation between the varying total magnetic fields 
on two sides of any sheet distribution of conducting material, irrespective 
of how these fields are produced. When there is a varying magnetic field 
of external origin, the total field will consist of this inducing field together 
with the field of the currents induced in the sheet. Let the potentials of 
these two parts be denoted by 2® and Q®, respectively, so that 

Q = 20400, 
Q© will satisfy Laplace’s equation and, in general, will tend to zero at 
great distances from the conducting sheet; Q© will satisfy Laplace's 
equation in the immediate neighbourhood of the sheet, and the corre. 
sponding field intensity (— grad Q) will be continuous on passing through 
the sheet. Since the normal component of the total field (—@éQ/én) is 
continuous at the sheet, it follows that the normal component of the 
induced field (—@Q®/én) is also continuous. Further 
H.—H grad(Q®Y —QO®), (10 

so that equation (9) may be written in the form 


p div grad(Q? —Q®)-+ grad p.grad(Q® —Q) 


@ (2Q0  aQo@ 
ot\ en on 


If the inducing field is known, 2 will be a given function of the time 





and space variables. To find the induced field it is necessary to determine 
the function Q® which satisfies Laplace’s equation, together with the 
appropriate conditions at infinity and the boundary condition (11) at the 
conducting sheet. When Q® has been found, the induced current distribu- 
tion may be obtained from (7) and (10). 


4. The non-uniform plane sheet 

When the conducting sheet lies in one plane, say the plane z = 0 in 
Cartesian coordinates, the symmetrical properties of the field of th 
induced current-system show that 


QO(x, y, 2) QO(x, y, —2). (12 
Hence, at z = 0, Q@ QM, (13 
and therefore Y = QW /27, (14 


so that equation (11) reduces to 


22Q)) ; 6 (EQ 2 (EQO' : 
“p nto + grad p.gradQ") +4 jet) a ~ Ie). 
a fy. ot e Fa at\ oz | 4 
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It is therefore necessary in this case to find the function Q® which 
i) satisfies Laplace’s equation in the region z > 0, (ii) tends to zero as 
- tends to -+-oo, and (iii) satisfies equation (15) at z = 0 if the sheet is 
of infinite extent. If it is of finite extent, (15) must be satisfied over that 
part of the plane z 0 occupied by the sheet and Q® must be zero or 
constant over the remaining part. 

When the sheet is uniform and of infinite extent, the second term on 
the left of (15) is zero, and the solution is very simple. We note here, for 
purposes of comparison, the solution when the inducing field is of the 
form Q(z-+-vt, x,y), where v is real or imaginary. When v is real we have 
the well-known case treated by Maxwell (5) of a magnetic system approach- 
ing the sheet with velocity v; when v is imaginary we have the case of a 
periodic inducing field, having a simple-harmonic time factor. Taking the 
inducing field as originating in the region z > 0, so that Q°-0 as 


-—o0, it is easily seen that the potential of the induced field in the 


region z > 0 is given by 

Q kQO(—z-+vt, x, y), (16) 
where k 2rv/(p+27v). (17) 
On the other side of the sheet (z < 0) we have, by (12), 

Q kQO(z+-vt, x,y), (18) 


so that the total field on this side is 





I~ 
t— 


(1—k)Q@ — —P ge, (19) 
p+2rv 

showing that the screening effect of the sheet is simply a reduction of the 

intensity of the magnetic field without any change of distribution. When 
is real (Maxwell’s case) the intensity is reduced in the ratio p/(p+27v). 

When v is imaginary (periodic field), the modulus of the complex number 

p(p+27v) gives the reduction in amplitude, and its argument gives the 

change of phase. 

These simple results cannot be extended to a non-uniform sheet, where 
the distribution as well as the intensity of the field will always be affected. 
To obtain an analytic solution in this case it is necessary to assume an 
ippropriate series or integral expression for Q®, say > J,(t)f,(v.y,z) or 

8 


8 


| L(t, u)f(w. y, 2, u) du, where the functions f, or f satisfy conditions (i) and 
li) above. The coefficients J,(t) or J(t,w) are then evaluated by substi- 
tuting in equation (15). This leads to a set of difference-differential 
equations to determine the coefficients J,(t), or to an integro-differential 
equation to determine J(t, w). 


Che case where a series expression is the appropriate form for Q0 is 
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treated in §§ 5 to 8 of this paper. An illustrative case which requires the 
integral form for Q®, and leads to an integral equation for J(t, ~), has beep 
discussed by A. A. Ashour in a paper now in preparation. 

A quite different method of obtaining the solution by successive numeri. 
cal approximations is briefly indicated in § 13. This method seems like); 
to be of considerable value in geophysical applications. | 


5. Plane sheet with p and 2 doubly-periodic functions of x and ; 

If p and Q© are periodic functions with the same period-rectangles, ther 
Q© will also be periodic. We take the case where 

p = 2mo{l—e(cos px+cosqy)}, with «<1 (20 

and assume that the inducing field is an even function of x and y, and origi 

nates in the region z > 0, so that its potential is expressible in the for 


co co 


e) y aos COS xf 221 9292)) 
(yo) = a En. COS MpxX COS NGY EXP{z,/(m*p*?-+-n2q")}, (21 


m=0 n= 
where the coefficients Z,,, ,, are given functions of the time; the coefficient 
E,o, corresponding to a constant term in 2, may be taken to be zero. 
The potential of the induced field for the region z > 0 will then be of 
the form 
QO Din.» COS mpx cos ngy exp{—z,/(m?p?-+-nq?)}, (22 


myn 


Ms 


iMs 


m=0n=0 


where the coefficients J, ,, are unknown functions of the time. This expres- 


mn 
sion for Q® satisfies conditions (i) and (ii) of the previous section. Substi- 
tuting the above expressions into the boundary condition (15), we obtain 


io @ x 
> > cel, (mp? sin px sin mpx cos nqgy-+-ngq? sin qy cos mpx sin nqy)+ 
m=0 n=0 
tol, ,(m*p?+ n?q")(1—e cos pa—e cos qy)cos mpx cos nqy+- 
| 
1 PAmn 2n2_1 m2n8\ang ial 
a i. (m*p*+n*q")cos mpx cos nqy 
: dE 9s 9.6 bp 
= a/(m*p?+-n?q")cos mpax cos nqy. (23 
tf 
m=0 n=0 


The expression on the left of (23) is easily reduced, by means of trigono- 
metrical identities and a rearrangement of the terms, to a double Fourier 
series, so that (23) may be written in the form 
co ice) 
© ¢ 6 © ] € € ' 9 9 
> cos mpx cos ngy[{o(m?p?+ n?q?) + 4/(m2p?+-n2q?)( d/dt) nn— 
m=0n=0 r 


Te i il 
—1oe{m(m—1)p?+n?q3F,,_1,—foe{m(m+1)p? 4+ Qn an— 
—oefm?p?+n(n—1)q7}0,, 1 —doe{mp?+-n(n+ 1)q? Inn val 
oe f+ a) . 
= > > cosmpxcos nqy,(m*p?+-n7q?)(d/dt)E, ns 


m=0 n=0 


(24 
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where, for the sake of generality in the notation, we introduce the two 
additional sets of coefficients J_,,, and J,,_,, all these coefficients being 
defined to be zero; this gives correct values for the terms corresponding 
to m 0 and to n Q. 

Since equation (24) holds over the whole period-rectangle, we may 
equate coefficients of cos mpx cos nqy on the two sides, thus obtaining the 
doubly-infinite system of simultaneous linear differential equations of the 


form 
{ 
d 
- he ln Cmyn si ie i+ ann mn din Riis +1 "m,n Reean 
dt 
l d 
: * = on 
> 2a Geom” (m i, © > (25) 
a.) (m*p"+-n-q )a 
to determine the unknown functions J,,,,. The values of the coefficients 
9s Om. noee IN (25) are readily found by comparison with equation (24). 


Since the double Fourier series on the right of (22) when z = 0 represents 
: function which is finite and continuous for all 2 and y, it follows that 
I, tends to zero when m tends to infinity and also when n tends to 
nfinity. These conditions together with the conditions that L,,,, J, -1, 
ind J), are all zero, are just sufficient to determine the required functions 
|, uniquely. The labour of solving the difference-equations (25) subject 
to the above boundary conditions is evidently considerable. Thus to 
obtain the Nth-convergents of the J,,,,’s (that is, the values of J[,,,, for 
m< N and n < N, on the assumption that [,,,, = 0 for all m > N or 
NV), it is necessary to solve (V+ 1)?—1 simultaneous equations. It is, 
however, possible to devise iterative methods in which the Nth-convergents 
an be obtained from the (V—1)th-convergents. This is illustrated by the 
relatively simple problem which we now consider. 


6. Plane sheet with p a periodic function of x 
Some important features of the induced current systems which flow in 

non-uniform sheets can be illustrated by considering a case where the 
inducing field is such that, in a uniform sheet, it would induce currents 
flowing in a direction which is perpendicular to the ridges of high resistance 
in the non-uniform sheet. One obvious feature would be a deflexion of the 
urrent lines by the ridges of high resistance, but other interesting features 
if the current system also emerge. We therefore consider the case where 
the potential of the inducing field is 

(J E(t)e% COs Gy, (26) 
ind the resistance of the sheet is given by 

p 27a(1—ecos px). (27) 

u 
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In this case we can evidently assume that the potential of the induced 
field is of the form 


ae Q - ~ 2 I 22) R 
QO = cosqy 3 m(t)cos mpx exp{—z,/(q?-+-m?p*)}. (28 


Following the general method indicated in § 5, it is found that the 
coefficients J,,(t) of the Fourier series on the right of (28) satisfy the 


m 
following recurrence relations 


I, = 2e-(1+ So (29 
(s iy =€ nig +(s+1)*D}L,—h, (30)t 
i Me! sco (nn*e- D fy In—In-1 (Mm > 1), (31 


m(m—1)s+1 ma = m—1)s-+-1 sit Ter oe 


ve 1 d 

where s = p?/q? and D denotes the operator — ae 
Between these equations we may sllehinene | 7 in succession to 
obtain a differential equation of the mth degree connecting J,,, J, and the 


known function H(t). Writing the equations in the abbreviated forms 


I, = a,e1),—), E, (32 
I, = 4, €712,,_,—),,f,,-. for m> 1, (33 


where the a’s and b’s are easily determined by comparison with (29), (30 
and (31), we find that 


} - Im(D)h—Pp(D)E, (34 
where Pn = G,, <*p,,_,+-5,, Pm-2> a by, Pe = 0, (35 
Ym = Im € = 1 +6,, Ym-2> = A *, VW = 1. (36 


It follows that p,, 
operator D, so that (34) is a differential equation of degree m. The ratio 
Pm/Ym iS easily seen to be the mth convergent of the infinite continued 


and q,, are polynomials of degree m in the differential 


fraction in the expression 


b, by bs 
—3 








tat 


0o— 2. 
\a, e1— a,e— age 


mS (37 


which may be obtained formally from (32) and (33) by writing (32) in the 


form ‘ b, E 
_ ay € a —T,/Iy’ 
and making successive substitutions for J,/J), 12/4... using (33). 
When z 0, the Fourier series on the right of (28) represents a function 
which is finite and continuous for all x. Hence J, > 0 as m>©. li 


m 


+ Note that (30) differs slightly from (31) with m 
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follows from (34) that J, must satisfy the limiting form of the differential 
oe Yn(D Ig = Py(D)E, (38) 


as m —> ©. 


7, Numerical solution for a periodic inducing field 

To determine the steady state solution when the inducing field has a 
simple harmonic time variation, we may take E(t) as the real part of 
£, e™ and replace the operator D in the above by ia/og. The value of the 
complex ratio J,/#, is then found from the infinite continued fraction (37), 
and the values of J,/#,, I,/F,,..., are determined in succession from (32) and 
(33). The potential of the induced field is then the real part of the resulting 
expression (28) when J,,(¢) is replaced by J, e*™. 

To illustrate with a numerical example, we take 


p q, sothat s=1, ) 


‘ (39) 
€ 3> xX oq — i, j 
so that the resistance of the sheet is given by 
Oma - 
p = —(l—§cos pz). (40) 
P 


With these values, the successive convergents to the continued fraction 
(37) for J,/E are 


Pi __ (-10000-+0-30000% ln 0-12853+-0-321122, 
4 de 
P3 __ (.13961+0-32171i Pa _ ().13315-+-.0-32160i, 
43 V4 
and Ps — 0-13321+0-32155i. 
Ys 


These are evidently tending to a limit which is not very different from 
P3/45- Taking this as the value of J,/E, the values of [,,/E, were calculated 


form = 0 to m = 4, and are given in Table | 


TABLE 1 





TIlE, Bm mod(Z5y/4 Ey ) | €m = arg(Iy/E,) 
13321 [55 0°345 | 67°5 
7908 7 O01 ph 514 
125 7 0°027 40°4 
1d 77 saa 31°4 
181 o°001g 14°9 





We deduce that a magnetic field, whose potential is 


Cyl E, e?* cos at cos py, (41) 
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will induce a current-system in the sheet, whose magnetic field in the 
region z > 0 has the potential 


2) 
Q = E,cospy S B,»cos(at+e,,)cos mpxexp{—pz,/(m2+-1)}, (4 la 
m=0 


where the values of f,, and e¢,, are given in Table 1. 
This result may be compared with that which would be obtained if the 
sheet were uniform with the mean conductivity 
27/p 
p da 3 
—=—s3 Dovel © ene ~~ ee” (43 
2x J 2mo(l—ecospxr) 2zovd 
0 


In this case the potential of the induced field for the region z > 0 is given 
by QO = 0-4082E, cos(at+ 66° )e-”* cos py. (44 

Since the stream-line function for the induced current is given by (14 
the current lines will be the curves Q = C, a constant, where Q( js 
obtained from (42) by putting s = 0. Also from (7) the induced current {| 
density is given by (a 
eee ee (45 

2a oy 2a Ox 

The current lines for the three epochs at = 0, 60°, and 120°, at = 360 
corresponding to a complete period of the inducing field, are shown in 
Fig. 1. At at = 180° the current lines will be the same as at at = 0, except 
that they are reversed in direction. During the second half of the period 
the current distribution will pass through the same sequence of changes 
as in the first half, but the currents will flow in the opposite direction. 
The distribution of currents, which would be induced by the same field 
in a uniform sheet of the same mean conductivity, are shown on the | 
right-hand side of the diagrams in Fig. 1. All the current lines are drawn 





for values of C increasing from zero by increments of 0-05. 

One obvious result of the non-uniform distribution of conductivity is 
the deflexion of the currents away from the region of high resistance, 80 
that they form circuits whose centres are, in general, in the region of low 
resistance. Fig. 1 shows that this is true at any rate for the greater part 
of the time. 

The current density reaches a maximum value when at = 114”, the dis- 
tribution being then practically the same as that shown in Fig. 1c. The 
minimum current density occurs when of = 24°, which is between the 
two epochs corresponding to Figs. 1a and 15; it will be observed that 
the current flow has changed in direction during the interval between . 
these two epochs. In the case of the uniform sheet the induced currents les 
simply decrease to zero intensity at time af = 24°, and then start to flow 
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again in the opposite direction. But, in the case of the non-uniform sheet, 
some remarkable changes of distribution occur during the period when 
the current flow is being reversed. To illustrate these changes, the current 
distribution was found for the epochs af = 17°, 24°, 30°, and 40°, and the 
corresponding current lines are shown in Fig. 2. These current lines are 
drawn for equal increments of 0-01 in C, except for the broken lines in 
Fig. 2a, which correspond to the much smaller increments of 0-0002 in ( 

Fig. 2a shows that the vortex-system of current, which is centred at 4 
on the line of highest conductivity, and which was of strength C,, = 0-240 
at time t = 0, has decreased in strength to C,, = 0-100 at time at = 17° 
Also a new vortex-system, of opposite polarity from that at A, has just 
come into existence at B, on the line of lowest conductivity. The strength 
of this new vortex is only —0-0004 at time at = 17°. 

At epoch at = 24° (Fig. 2b) the vortex-system centred at A has decreased 
still further in strength to the value 0-040. It has also shrunk in dimensions, 
its outer boundary being a rectangle which is gradually closing up in the 
x-direction. The new vortex-system, which came into existence at B at 
time at = 17°, is now seen to be actually a double system with a pair of 
centres which have now moved apart. The strength of this vortex-system 
has increased (numerically) to the value —0-032. Evidently there will not 
be in this case a particular instant at which the current density is every- 
where zero (as there is when the sheet is uniform), because, while the 
vortex at A is decreasing, the one which started at B is increasing. The 
exact epoch of minimum total current occurs very shortly after at = 24, 
at the instant when the strengths of the two vortices (of opposite polarity) 
become numerically equal, their values being in fact about + 0-036. 





At epoch at = 30° (Fig. 2c) the vortex-system which was centred at 4 
has closed up and disappeared completely, while the new twin vortex- 
system has increased in strength to —0-057, and has spread out so as to 
occupy the entire region. The two.centres of this new system have also 
moved still farther apart. Most of the current still circulates right round 
the pair of centres, which are of the same polarity, but some of it now 
passes through the region of high conductivity (cf. the current line 
C = —0-01), to link up with the adjacent vortex-systems. 

At epoch at = 40° (Fig. 2d) the greater part of the current is flowing 
nearly parallel to the 2-direction. But the most noticeable feature is now 
the complete separation of the twin vortex-systems which originated at B; 
one of the vortex centres is approaching the point A, and the corresponding 
vortex-system is beginning +o coalesce with an exactly similar system 00 
the other side of A. The centres of these two systems will continue to 
move closer together and will eventually coincide at A. The current 
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distribution will then be as in Fig. 16 (epoch at = 60°) where the centres 
have nearly but not quite coincided. 

The changes which take place in the current distribution during the 
general reversal of current flow can now be summarized as follows. Firstly. 
new vortex-systems of current, of polarity opposite from the old ones, 
appear in the regions of low conductivity. These new systems each split 
in two, and the new vortex centres move apart and migrate towards the 
regions of high conductivity. Pairs of vortex-systems coming from adja- 
cent regions of low conductivity approach each other and coalesce in the 
regions of high conductivity. Thus new vortex-systems are built up in 
the regions of high conductivity in which the currents are flowing in the 
direction opposite from that of the previous systems. 

It is evident that similar changes of current distribution will occur in 
any non-uniform sheet if the inducing field is periodic or fluctuating, since, 
when the field changes from an increasing one to a decreasing one (or vice 
versa) the new currents which are induced will be opposed by the old 
currents which linger in the regions of high conductivity. Consequently 
the new currents will be deflected so as to form closed circuits, which are 
initially in the regions of low conductivity. 

The time taken up by the above changes in the current distribution is 
about one-eighth of the half-period for the particular case considered 
above. In general this time will depend partly on the difference between 
the maximum and minimum values of the conductivity, and partly on 
the rate of decay of free current-systems in the sheet, which in turn is 
largely determined by the mean conductivity. The decay of free current- 
systems is considered in the next section. 


8. Aperiodic inducing fields: the effect of sudden changes 

If the inducing field changes instantaneously by some finite amount, 
and afterwards remains constant, there will be induced in the sheet a 
system of currents, which subsequently decays freely. To examine this 
vase, we take the time-factor H(t) in the expression (26) for the inducing 
field to be of the form E(t) = EH(t), (46) 


where £ is a constant and H(t) is Heaviside’s discontinuous function, 
defined by H(t) = 0 when t < 0 and H(t) = 1 when t > 0. The field of 
the induced currents is then given by (28), where the time-factor J,(t) is 
given, from (38), by 
I(t) = Elim 2™) AW), (47) 
m+ Ym(D) 


which can be evaluated by the methods of the operational calculus. The 
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maining time-factors /,,(¢) in (28) may then be found from (34) with the 
ist term on the right of this equation taken as zero. This process is in 
fact equivalent to solving the differential equations (29) to (31) for J,,(t) 
E and (0) = 0 for m > 0. The result 


m 


with the initial conditions J,(0) 
therefore corresponds, as it should do, to the free decay of currents in the 
sheet, the initial distribution of these currents having a magnetic field 


given DY QQ Ee-® cos qy. (48) 


When the effect of a sudden change in the inducing field has been found 
n the way indicated above, the result for any other aperiodic variation 
f the field can be deduced by the application of well-known theorems in 
the operational calculus. 

To correspond with the case of the periodic field already treated in § 7, 
calculations have been made for the aperiodic case when the distribution 
f conductivity in the sheet is given by (40), and with p = q. To determine 


L(t), the successive convergents 





»P,(D 
I,,(t) = BP) AW), (49) 
; qn(D) 
where nm = 1, 2...., to the expression (47) for J,(t) were calculated, using 


Heaviside’s partial fraction rule, since p,(D) and q,(D) are polynomials 
each of degree » in the operator D = (1/oq)d/dt. The first three con- 
vergents are respectively I 


oa = Be-om, (50) 
I, = E{0-7316e-075990a 4 ()-2684¢-1'65430a0) (51) 
fio E§0-5977¢ 0°69460qt _| 0-387 Lle—141200gt | ()-)1 52e—2'544ea | (52) 


Although the above expressions for I, and I,, are quite different, it is 
remarkable how little the values of I), and J, differ throughout the range 
'=(0 to t = 4/oq, by which time J,(t) has decayed to 0-039 of its initial 
It is 
evident that J, will differ still less from J) 5; we therefore assume that J), 


value. This will be seen from the second and third lines of Table 2. 


is a good approximation to J,, and we use (52) together with (34) to 
calculate J, and J,; to this degree of approximation, the values of J,, for 


m > 2 are all zero. The values of J, and J, are shown in the last two lines 
f Table 2. 


TABLE 2 








O°2 O*4 06 0°38 Io | 2°0 3°0 4°0 
in oe a > ee | j o 
2£|| 0-9054 8212 | 06784 5632 | 0°4699 | 0°3935 | 0°1698 | 0°0767 | 0°0350 
( ) 8209 | 0°6783 5633 | 0°4698 | 0°3939 | O1721 0:0800 | 0:0386 
} 741 | o-1I74 1406 | O°1503 | O°1513 | O°1077 | 0°0612 | 0°0324 
2 00 44 | 0:0008 0158 | 00207 | 0:0241 | 0°0248 | 0-0162 | 00090 
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Med geees 


(7 = 


time 
from 


We conclude that a sudden change at time ¢ = 0 in the inducing field 
of amount represented by the term 


in its potential, will give rise to a system of currents in the sheet, whos 
potential in the region z > 0 is given, for any subsequent time, by 


where J,, J,, and J, are given in Table 2. 

The current lines at any particular instant are given by QY = C, wher 
Q® is obtained from (54) on putting z = 0. Fig. 3 shows the current lines 
corresponding to equal increments of 0-1 in C, for the three epochs 
ogt = 0-2, 1, and 2. This figure shows the expected deflexion of th 
currents by the ridges of high resistance {x = (2n+-1)z} where n = 0, | 


(~ = 2n7) of maximum conductivity. 

The total current circulating in the sheet at any time is determined by 
the maximum value of Q” at that time. From (54) the maximum valu 
of QY at time ¢ is J,(t)+4,(t)+1,(t), and is attained at the vortex centres 


long period, the time taken up by the changes in distribution, described 
in § 7, which occur during the general reversal of flow of the current, will 
be roughly of the same order. Thus, in the case worked out in § 7, the 


9. Thin spherical shell with non-uniform conductivity 

We now consider the case of a thin spherical shell of radius a for which 
the resistance p at any point is a function of the angular coordinates (6,A 
The equation (11) in polar coordinates (r, #,A) becomes 


io @) n 
where { = r/a, > denotes the double summation } } and £7, Im are 
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Q© = HH(tle“cosqy (z > 0) (53 


QO = {he-’+I, e-®*? cos qx+ I, e-®** cos 2¢x}cos qy, (54 32 


so that they form closed circuits centred at points on the lines 


2u7, y = 2n7). Hence from Table 2 we see that the total current 


will be reduced to about 0-3 of its initial value in the time t = 2/oq. We 
may conclude that, in the case of a periodic inducing field of sufficient] 


taken was about one-sixteenth of the complete period 2z7/«, and 
(39), « = oq/3, so that the time taken is approximately 1-2/0q. 





l (0p @ 1 dp a 
asatan na 
a?\ 00 20 ' sin?@ dA eA 





P o ai 92 p # Vow _ow 
+365 (sn a) t aoe al oe 


: : 


= — An — (Q040%) at r=a. (55 
oter 


Since Q® satisfies Laplace’s equation in the neighbourhood of the shell 
we may analyse it into spherical harmonics in the form 


QO =a > (¢” Em+[-" -1]m)eimA Pm (cos 6), (56 


m,n 
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Fic. 3. The freely decaying current ‘system produced initially by a canine, Ce ~~ = 
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complex numbers, their real and imaginary parts being known functions 
of the time. The real part of the expression (56) can be made, by suitable 
choice of Ey and J”, to represent the potential of any given inducing field, 
The E£% terms correspond to the part of the inducing field arising from 
sources external to the sphere of radius 7, while the /%” terms correspond 
to the part of internal origin. We shall assume that the field originates 
outside the shell so that near ry =a the J™’s will all be identically 
zero. . 

In the same way, since the potential Q® of the induced magnetic field 
satisfies Laplace’s equation, and originates from currents in the shell. it 
can be expressed in the forms 


QW —a y ime-n lgimA P™(cos 8) (r - a), (57 
m,n 

Ow —a > emeneimA P™( cos 8) (r < a). (58) 
m,n 


The coefficients 77” and e?, which occur in these expressions, are complex 
functions of the time, and the problem is to determine them in terms of 
the known time-functions 2”. 

Since the normal component (—@éQ/ér) of the induced field is con- 
tinuous at r = a, we have that 

> (n+1)imeimP™(cos 0) = — Y nemeim Pm(cos 8), (59) 
m,n mn 
and, since this is true over the whole sphere, we may equate the coefficients 
of the separate harmonics, giving 
(n+1)™ = —mne™. (60) 
Hence, from (57) and (60) we have that 
O® = —a > ; 1" _ emg—n-lgim\ Pm eos 6) for r>da. (61) 
an N+ 1 
m,n 
which gives the induced field outside the shell in terms of the induced 
field (58) inside the shell. 

Substituting from equations (56), (58), and (61) into (55) we obtain 

wet ee wasn) + 
00 26 ' sin?@ @\ @A * sin 8 20 O} 














= & | 2n-+-1 MpimA Pm 36 
1 sim ar2| p3 n+l en‘ n (COS 6) 
dk™ dem ; 
— In nea | n \ pimA Pm Q 62) 
P z n( dt dt Ee wyeer). © 
mn 
This equation may be simplified by —_ the relation 
] 1 as 26 
1 — 0), 63) 
sin 6 als Ot a ae vilaiine \ 
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which is satisfied by any surface harmonic of degree n. The equation (62) 


is thereby reduced to the form 


nh ' 
- 
mn € 


ie —=/(d, d ; . 
SY eV se mA P™/(cos 0)! 47a > (7 Be + 562 }ne™ Paicos), (64) 
mn 


where V, is the operator defined by 


a 2n-+-l{ép ¢ 1 ¢ép 


5 
: -|-—— ——n(n-+1 (65) 
n+1 106 00° sin?6 0A 00 Ip 


f 
This operator depends on the distribution of conductivity in the shell; it 
reduces to a constant when the distribution is uniform. 

The equation (64) is satisfied over the whole surface of the sphere. 
Hence, multiplying by e-’“4P¥(cos@), and integrating over the surface 
f the sphere, we obtain 


4 em | e-iMA PM (cos 8)V, {e™4P™(cos 8)}sin 6 dbdA 
. 0 0 


i 


tra S 7 | “ Em .. n | | | elm—M)A P™( cos 0) PY (cos 8)sin 6 dédA 


hot dt 
. 0 0 
a a d (N+M)! 47 
4rraN E™M 4 eM : < 66 
& VT aH | (N—M)! 2N+1 rn 


by the orthogonal properties of tesseral and zonal harmonics. 

The coefficients of e” in the expansion on the left of equation (66) are 
constants which can be determined if the function p is given. Hence 
equation (66) is a linear differential equation in ¢ with constant coefficients. 
There will be an equation of this form corresponding to each pair of 
positive integers (M, NV) where VN > M > 0. The required time functions 
e” are therefore determined by a doubly infinite set of simultaneous linear 
differential equations (66). 

It will be seen that, if p is a function of @ only, the coefficient of e7 in the 
expansion on the left of equation (66) will be zero unless m M. Hence 
the functions e” which have different values of m will be independent ; 
this corresponds to the fact that the harmonics of different rank in the 
potential of the induced magnetic field will be independent. In this case 
we obtain a singly infinite system of equations to determine each set of 
coefficients corresponding to a given value of m. 

The solution of these equations will now be considered for a particular 
case of the function p. 
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10. Spherical shell with resistance p = p,(1+ cos @) 
Taking the simple distribution of resistance 





Pp = po(l+ecos 8), (67 
the operator V,, becomes 
2n-+-1 ss, ' 
Vv, =— = 7 Poeen —— n(2n-+-1)po(1+-e cos @). (68 


Hence, using the relations 
= 
(2n+ 1)sin =, PR = n(n—m+1)PM,,—(n+1)\(m+n)P™,, (69 
(2n+-1)cos@ P= (n—m+1)P™,,+(m+n)Pr,. (70 


which are satisfied by the associated Legendre functions, we find that 


nm+2 
Im __ > ») 
V, PS = ar po en(n—m-+- 1) P™, ,—(2n+-1)np, P™ 
—(n—1)(n+m)p,¢P™,. (71 
Instead of using the general result (66), it is simpler in this case t 
substitute directly into equation (64). This gives 
9 
ima} 21 4 > 
-> eme™ jar Ty Pe en(n—m+1)P™, + 


m,n 
+(2n-+-1)np, P™+-(n—1)(n-+ m)p,eP™ | 
} 


; 


= 4ra z (5 Em +5) neimPm (72 


m,n 
Since each side of this equation is an expansion in tesseral (or zonal 
harmonics of rank m, and the equation is satisfied over the whole surface 


of the sphere, we may ecuate the coefficients of harmonics of the same 


degree. Thus, picking out the coefficient of e’”AP™, we find 


Po e(n?— 1)(n—m)e”_, + {ro n?(2n-+-1)+ tran’ | em 
t 


d ao 


2 2? p 
+ py en*(n+m-+1)em,, = - 4nan*— E™, = (i 


Corresponding to each integer m > 0, there is an infinite set of these 
equations extending from » = m to n= oo. It will be seen that the 
correct equation is obtained in the special case n = m, since the first 
term on the left of the equation then becomes zero. 

[f the inducing field contains only the single term E”, the above equa 
tions can be solved by methods exactly similar to those used in §§ 6 to $ 
the expression for e” reducing to a continued fraction similar to (37). 
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In the general case when all the #”’s are present, the equations begin- 


ning at m may be written in the form 


m lpm . 1m 
m+ Ag € “emt koe ER: 
’ lpm - -1 7 7 = 4\ 
em+s Og-y € CM ga t+ Osa Cmte—-at hee Eis (8 > 1), (74) 
shere the coefficients a), ky, @,_1,--- are linear functions of the operator 


) =d/dt, which are easily determined by comparison with (73). By 
elimination of the intermediate coefficients from (74) it is possible to 


express e”.. in terms of e” and the known functions E™ ,, for r = 0 to s—1. 


m “m 
Writing 8-1 
m ae me _ him yA 
em+s Gs —m > Prs kr 4m+r? (75) 
r=0 


ve find that g, and p,, satisfy the recurrence relations 
1 "a 
TV a. 1€ Vs-1 b. 1 Vs—2> qi Ay € > qo s (76) 
1 — — 
D.. Ae—1 © Pres b. 1 Pr,s—2: Prrti a Prr = 0. (77) 


The coefficients g, and p,, may also be expressed as continuants in the form 


Ys Ko, 1\€); Pr,s Ky +1,s-1(€): (78) 
where 
ae} | 0 0 
1 
b, 1 a, 1€ ] 0 
K.(e) > thn aue*’-. . i (79) 
1 
i a,_1€ —] 
0 0 0 : b. a,<«-* 


This shows that q, is a polynomial of degree s in the operator D, and p,., 
is a polynomial of degree s—r—1 in D, so that (75) is a linear differential 
equation of order s connecting e™,, and e™. 

Since the expression (58) for Q® represents a function which is finite 
and continuous over the sphere r a, it follows that e™,,—> 0 as s > 0. 
Hence from (75), e” must satisfy the limiting form of the differential 
equation — 
q,(D)e; > p,(D)k, ebm, 0 (80) 


7 0 


as soo. We therefore write 


s—1 
em lim— S a (81) 
s onl — 
r=0 


and obtain successive convergents to e” by using (76) and (77). When 
the field is periodic with period 27/x« we may replace D by ia and follow 


the same procedure as in § 7. When the field is aperiodic, the successive 
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convergents to (81) may be evaluated by using Heaviside’s partial fractioy 
rule, in a way similar to that adopted in § 8. 


When e” has been found, the remaining coefficients e”,,, €”..,..., mal 


m+1> “m+2 
be obtained in succession by using equations (74). 
Some numerical calculations of this solution for the particular case whe; 


the inducing field is of the form 
Q© = al (t)cosAsin 0 (89 


have been made by A. A. Ashour and the writer (4, 1948) and have bee 
used to discuss some effects of the non-uniform distribution of condue. 
tivity in the ionosphere. 


11. Solution in ascending powers of « 

Instead of finding successive convergents to the limiting expression (81 
for e” in the manner indicated above, it is also possible to express ¢ as 
a power series in ascending powers of e?. To do this we first obtain th 
‘corresponding’ series for the limiting ratio K,,,(e)/Ko,.(€) as s > 00, wher 
K,.(e) is the continuant defined in (79). 

Expanding K, (ec) by Sylvester’s rule, we find that 


K, {€) = e*-*—1a,4,,,...4,(1+- > A,,€*), (83 
where A,,=0 for 2n>s—r+1, (84 
S—2n+2 s—2n+4 s 
and A,, > > - D>  B,B,---B, for 2n<s—r+1, (85 
h=Pt1 t=ht+2  tn=t—it2 
where B, = 6,/a,_, %- (86 
K r kk. J 2n 
Hence, rasl€) it seems - é Aire (87 
° ” b ] 
Ky(€) py 4y...0,_1 1+ > Ange 


and provided all the a’s and b’s have finite non-zero values, this ratio is 
an analytic function of « (regarded as a complex variable) near the origin 
We may therefore expand it by Taylor’s theorem in the form 


q r ve ‘ 
K,,e(€) — : ( £3 s Cc e2n ), (SS 
K,,(€)  @@;...2,_, — 
where (1-4 > c,e2")(1+ p A,,,e™) = ]+ > A alee (59 


Equating coefficients of like powers of e? on the two sides of the identity 
(89), we have that 


A (90) 


Cptlpy Ajo tly Agqt+--+Apo Br 


Dp 


. ie : m : : 7 vo And 
for all positive integral values of p. Taking p = 1 and using (85) we find 
that 


s 
Cy A,,—Aj = a. B, — 2 B, 
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| fractin Taking p 2 we have similarly that 
Cy A., hee C,Aj9 
dyeees I e—_? g e—9? ¢ r s 
s > BB ; = . = > 
ya fe PP > a B,, Br, } 3 “th > B,, 
r+1t:=h+2 th=1t2=+2 =1  t3=1 
ASE Whe 
4 . < >< 
() QO Q”) 
how ya Mt; Mts hn ht Bi, Bre hn ht B,, B,,: (92) 
1=1t=h+2 ti=1t:=1 t=1t:=1 
It is now easy to prove by induction that, for all integral values of 
ave be l(s—r+). 
cond j r &+1 K-93 " 
c (—1)"> >... > &,B,---f- (93) 
M41 090 t= 1 
For values of n }(s—r-—1) the coefficients c, will depend on s, but 
3) shows that for n }(s—r-+1) the coefficients c,, are independent of s. 
S10n (Ss ° . . ' . 
Hence, if we make s > o the series on the right of (88) defines a power 
series which has the property that its terms as far as e?"*" for any positive 
tain ; ca) : os , — ‘ 
nteger n, are identical with those of the Taylor expansion of K, ,(e)/ Ko ,(e) 
whe , 
lor every s 2n+r—l1. 
In the special case r |. the ratio A, ,(e)/Ko,(e) is the (s+1)th con- 
vergent of the infinite continued fraction 
| h b, e be hoe 
| 1 2 at L 2 (94) 
R4 A, € 1 ay€ l ly € 1_) Xo | ay ! As 4 
1 (85 ind the existence of a series having the above property is well known 
Perron, 6, 1912). It is called the ‘corresponding’ series of the fraction, 
86 ind the relation between the fraction and the series is indicated by 
e b,€* db, €* € — 9 F 
) <——_...~ (1 > C,¢" 1. (95) 
Ag+ 4+ a2 lo n=1 
ratio is \dopting the same notation for the general case (i.e. for any positive 
> Origil ntegral value of 7), we write 
‘ K (é)} é’ 4 . 
lim — ~ {1 i. > c,e"1, (96) 
96 P Ke (€) Ay Ay---A,_4 n=1 
where ¢c_ is now given by the formula (93) for all values of n. 
89 \s far as the writer knows, the explicit formula (93) for the coefficient c,, 
dealt ff the corresponding series. either for the special case of the continued 


Iraction or for the more general case of the ratio of two continuants, has 


1 


ot previously been obtained, although, in the case of the continued 
lraction, both Stieltjes (7, 1889) and E. T. Whittaker (8, 1915) have 
we fl indicated general processes whereby any particular term may be calculated. 
Thus Whittaker has shown in effect that the coefficient C. 


eading element of the 2nth power of a certain matrix. 


9092.7 


is equal to the 
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Xx 
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The sign ~ cannot in general be replaced by the sign of equality, sinc 
the process of deriving the ‘corresponding’ series involves changing the 
order of a double limit, but if the expression on the left of (96) converge 
uniformly to the limit L(e) in a closed domain PD containing the origi 
then the series on the right of (96) is convergent inside any circle ¢ =, 
within D, and its sum is L(e). This may be proved by a simple extensio) 
of the corresponding theorem given by Perron (6, chap. 8, theorem 20) fo 
infinite continued fractions. The uniform convergence of the expressio; 
on the left of (96) follows in the present case from the fact that a, and! 
tend to finite limits (namely —2 and —1) as s > ~, and we are therefon 
justified in equating the two sides of (96). 

We can now use (96) together with (78) to express the limiting value o} 
P,s/Us AS 8—> oO in the form of a power series in €*, explicit formulae fo 
the coefficients of this power series being obtained from (93). On subst 
tuting in (81) we get a corresponding expression for e”” as a power serie: 
in e*. The coefficients of this power series will be rational functions of th 
operator D, and the evaluation of e” as an explicit function of the tim 
for both periodic and aperiodic fields can be effected by methods similar 
to those of the previous sections. 

In practice it has been found that the method of solution described ir 
§$ 10 is generally more convenient than the one just described, but the 
latter has the advantage that one calculation serves for all values of « 
the series obtained being fairly rapidly convergent if « < 1/3. 


12. Applicability of ‘relaxation’ methods to the problem 

The analytic methods described in the above §§ 4-11 for obtaining solu- 
tions of the fundamental equations of § 3 lead to manageable calculations 
only when the distribution of conductivity is of a fairly simple form. In 
some geophysical problems it is desirable to estimate the effects o! 
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irregular distributions of conductivity, given as empirical functions ove! 
the surface of the earth (for example the integrated conductivity of th 
oceans), and to deal with these cases other methods for approximating to 
the solution must be found. 

The possibility of using ‘relaxation’ methods for this purpose is firs! 
considered. The present problem involves the determination of a potentia 
function Q® which satisfies Laplace’s equation, with appropriate condi- 
tions at infinity, and the boundary condition expressed by equation (I! 
of § 3 at the conducting surface. The latter boundary condition involves 
the time-variable, but in the case of periodic fields (to which we shall now 
confine our discussion) this variable can be eliminated if (2 is treated as 
complex, as in § 7. The required function © is then, in general, a complex 
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function of three space-variables. The familiar relaxation methods for 
slving Laplace’s equation can be used to determine Q” provided the 
number of space-variables can be reduced to two, but this is only possible 
when the distribution of conductivity is expressible in terms of only one 
variable. For example, in the case of the spherical shell, if the conductivity 
isa function of 6 only, the inducing field Q® may be analysed into spherical 
harmonics of the form E”P”(cos@)cos(mA+<,,), and the harmonics of 
different rank (m) treated separately, as shown in § 10. The term in Q” 
corresponding to the harmonic of rank m will then be of the form 


Q"" cos(mA-+t-e,,), 


i 


where 2” is a (complex) function of the two variables (r, 9), which can be 


m 
found by relaxation methods. 

The above shows that the straightforward application of relaxation 
methods makes it possible to extend the solution obtained in § 10 (for 
which numerical calculations have been made by Ashour and Price) to 
the case where p is any function of @. Relaxation methods alone will not, 
however, enable us to deal with the case where p is a function of both @ 
and A, unless such methods are extended to three dimensions. This, of 
course, would be theoretically possible, but does not seem likely to lead 
to a convenient method for dealing with the present problem. What is 
really required is some way of analysing the mathematical problem into 
1 number of simpler ones, each of which involves a manageable set of 
calculations. Two possible ways of doing this, suggested by physical 


considerations, are described in the next section. 


13. Alternative methods of approximating to the solution 
If the equation (11) is written in the form 
; : QO aro = 
p div grad ‘t+ grad p.grad ‘¥ = ——__- — —_., (97) 
oton cton 
where ¥’ is the current function given by (8). it will be seen that ‘Y can 
be divided into two parts ‘ and , which respectively satisfy the 


equations 


A2( ye) 
. . _ Ose 
p div grad 4 grad p.grad Yo — ———., (98) 
cton 
; Para yr) 
p div grad ‘V+ grad p.grad YO — — —_., (99) 
oton 


The first part “ corresponds to the current system which would be 
obtained if self-induction were ignored. This part can be calculated 
directly from (98) since Q® is a given function, the calculation being made 
by relaxation or other numerical methods, if analytic methods are not 
convenient 
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The second part ‘ represents the effect of self-induction. It cannot 
of course, be calculated directly from (99) because Q® is unknown. By 
under certain circumstances (which are examined below) a first approxima- 
tion Q, to Q” can be obtained by calculating the magnetic potential of the 
current system ‘. In the case of the plane sheet or spherical shell this 
can be done by known methods. If we now substitute Q, for Q on the 
right of (97), and solve the resulting equation by ‘relaxation’ or other 
methods, we get a corresponding approximation ‘, to ¥’. 

A second approximation ‘Y, can now be obtained by calculating the 
potential Q, of the field due to 4, substituting Q, for Q® in (97), and 
again solving this equation. 

This process can evidently be continued and, if convergent, will lead to 
the required function ‘’ and the corresponding potential Q. In practice 
it is slightly more convenient to use (99) instead of (97) in the iteration 
process; this leads to the sequence of correction terms ‘Y, —‘Y, ¥,—%,.... 
which must, of course, tend to zero. 

To determine the conditions under which this process will lead to a 
solution, it is sufficient to consider a trivial case for which the analytic 
solution is known. We consider the case of a plane uniform sheet, so that 
equation (11) reduces to 


er ey 200 @Qo 
— + —_} = —-— incon: i 0. (100 
\0x* — — ay* otoz otcz 


and the relation between ‘’ and Q® is given by (14). If the inducing field 


is given by QO = ecos a cial (101) 


the analytic solution for the induced field (ef. (16) and (17)) is given by 


oo — _B 2 rz cos Aa ei ™, (102) 
1+B 
where B = 2zia/Ap. (105 
Now applying the method described above, we find that 


, bow , 
yo — T, cosa el, 


\p 
whence (Q,).<9 = 27'f = Boosdaxe™, (104) 
and therefore Q, = Be-*“ cos Ax ei, (105 


In the next step the equations are of exactly the same form as before. 
except that Q, replaces Q; it easily follows that 


Q,—Q, = —fe-** cos Ax et, 
and continuing the process, we evidently get 


Q,, = B{1—B+f?—...+(—B)"-NHe- cos Ax ef. (106) 


n 
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Hence this process converges to the correct solution (102), provided that 
(107) 


) 


It will be observed that 8 depends on the conductivity (1/p) of the sheet, 
and on the period (27/«) and spatial distribution (wave-length 27/A) of 
the inducing field. 

The above indicates that this method will in general work best for the 
case of low conductivities, slow oscillations, and high (spatial) harmonics 
in the total field. The influence of the spatial distribution of the field is 
particularly important in the case of non-uniform distributions, for reasons 
indicated below. 

The method just described is based on the possibility of obtaining a 
satisfactory first approximation by neglecting self-induction. This method 
is not satisfactory if the conductivity is too high, but physical considera- 
tions suggest that in this case a suitable first approximation may be 


obtained by neglecting the resistance. Thus, putting p = 0 in (97), and 


denoting the resulting approximation to Q® by Q® we obtain 
EQ EQ 
- (108) 
on cn 


from which 2” may be determined since Q® is known. In the cases of the 
plane sheet and spherical shell the formulae required for this determination 
are well known. 

Corresponding to this value of Q” a current function ‘Y can be found, 
and substituting ‘Y in (97) a second approximation Q® to Q™ is obtained. 
If this procedure is repeated, a third approximation 2 can be obtained, 
and so on 

Applying this method to the simple case already treated by the first 


method, we easily find 


n—1 
G2 I i inc \, Az eos Aa ei, (L109) 
| 5 5 B | 
which tends to the correct limit (102), provided 
Zara 
3 = © (110) 
Ap 


Hence comparing (110) with (107) we see that the two methods are 
complementary, in the sense that when one does not converge the other 
dicate con theialiens expect one or other of these methods to be 
ipplicable in all Cases. 

A point arises in connexion with the application of the second method 


to non-uniform distributions, which is not brought out by the above 
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example. This relates to the spatial distribution of the computed field Q” 
This distribution remains constant in the above example, but in the genera] 
case it will change at each step in the approximation. Now, the inequality 
(110) shows that for the second method to be convergent the field mus 
not contain spatial harmonics of excessively high order. For example, i 
the case of the plane sheet the wave-length 27/A must not be too small 
and, correspondingly, in the case of the spherical shell, the spherical 
harmonics in the field must not be too high. While this condition may ly 
fulfilled by the inducing field itself, and possibly by the first few approxi. 
mations Q®, Q@)...., to Q® there may come a stage when 0”) contains 
harmonics of such high order that the method is not convergent. If this 
occurs, however, it is always possible to separate the part of the field 
containing the high harmonics, and deal with this part by the first method 
Hence for some problems a combination of the two methods will be required 
[t will be evident that physical considerations can be of great assistance in 
guiding the work and reducing the labour of computing such problems. 
The two methods described above have been tested by J. M. de Wet 
(9, 1948) against the analytic solution for the particular non-uniforn 
spherical shell problem treated by Ashour and Price (4, 1948), and are 
found to work well when proper precautions, such as those noted above, 
are taken to ensure the convergence of the method. It seems likely that 
these methods will prove useful in considering the effects of the oceans 


on geomagnetic variations. 
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ON WEBER’S FUNCTION 
py C.G. DARWIN (National Physical Laboratory, Teddington, Middlesex) 


Received 2 November 1948] 


SUMMARY 


d=u 


The equation a3 bx®—a)u 0 occurs in a number of wave problems, and 
les of its real solutions are in course of preparation. The solutions are specializa- 
ns of the I function of Whittaker, but the present work proceeds from first 
rinciples so as to avoid this over-generalization. The appropriate solutions are 
or images 1n 2 0. For positive a the first of them resembles an exponential 
r the origin, rises to a large value near 2va, and then oscillates with slowly 


creasing amplitude and wave-length. The second resembles an exponential with 


} 


vative argument and after 2va also oscillates. Near infinity the two curves are 


\)° out of phase. For negative a there is a similar pair of solutions, but the oscilla- 
ns start at zero, and are roughly of equal amplitude ; near infinity they are 90 
if phase 
Convergent power series are given for the two solutions, and also the associated 
symptotic series for large x. Neither of these types of series would be much use 
nais large, whether positive or negative, and another set of series is developed 
plicable for large a. For positive a there are two such series, one applicable for 
ies of x from infinity down to somewhere not far above 2va, and the other for 
ies of a fi zero to a point somewhat below 2va. For large negative a the 
lar series appears to converge for all values of x. 


1. THE purpose of the present note is to discuss the various solutions of the 
lifferential equation d2u 
+ (fa?—a)u 0. (1.1) 
dx 

[he equation has occurred in various problems of vibration, and in par- 
ticular in connexion with the propagation of electro-magnetic waves in the 
onosphere. Tables are in course of preparation for its solutions; they are 
ing prepared by Scientific Computing Service Ltd., on behalf of the 
National Physical Laboratory, and it was this that stimulated the under- 
taking of the present work 

The equation is discussed at some length in Whittaker and Watson's 
Vodern Analysis (§ 16.5), but with the difference that the term in 2? is 
negative, signifying that the real axis is 45° different from that taken here. 
rhis, of course, makes no difference to its character in the complex plane, 
ut makes a very considerable difference when the question of real solutions 
with real a is considered; for example, the scales of relation, so valuable 
when x? has negative sign, are no longer serviceable. The approach in 
Whittaker and Watson is by specializing a more general function, and it 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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does not fully consider which among the various solutions are practical); 
most convenient. This approach from the general Whittaker function wa 
also used by J. C. P. Miller and C. W. Jones in the preliminary work fy; 
the forthcoming tables, but it of course entails mastering the properties 
of more general forms than those actually required. The present work has 
therefore started from first principles, though broadly following simila) 
lines. At the end, however, certain asymptotic solutions are develope 
which are not to be found in Whittaker and Watson. 

The accepted form of the equation contains the term }2?, but it is }y 
no means clear why the } was introduced. It complicates the algebra a 
all the solutions and disfigures nearly every formula by an unnecessar 
factor v2, which would not have occurred if the term had been take; 
simply as 2. In fact most of the present work was done without th 
factor }, but I have concluded with some reluctance to restore it so ast 
agree- with the form which is to be used in the tables, where, after du 
consideration, it was decided to retain it. If the function were to have th 
same sort of widespread use as the Bessel functions it would certainly by 
best to alter it, but with its more restricted importance, and with numerica 
tables which are to appear in this form, it seemed that the inconvenience 
of changing over would outweigh the inconvenience of the accepted forn 

The general geometrical character of the solutions is easily seen. First 
take the more interesting case of a > 0. From (1.1) it appears that 
x| < 2va, the curve is convex to the a-axis, while if |a! > 2va it 
concave. It is clear that there are two solutions, one an even function, an 
the other odd, and these make a convenient starting-point. Near th 
origin the even function, say wp. is like a hyperbolic cosine. It graduall 
drops below the curve coshava, but rises to a rather large value at 2) 
which is a point of inflexion. It then turns gradually towards the axis 
crosses it and executes oscillations, slowly diminishing in amplitude an 
wave-length. The odd function wu, starts like a hyperbolic sine and ther 
behaves in the same way, bending back and oscillating after 2va. Fo 
large positive values of x the two solutions become nearly proportional t 
one another, while for large negative values they are nearly proportiona 
but of opposite sign. This implies that these two solutions are not th 
best to take, but that instead one solution should be a combination whic 


cancels out the main asymptotic expression; in fact they should cor 


spond to the exponential functions instead of the hyperbolic. From thi 


form of the differential equation it is evidently possible to take the tw 
solutions as mirror images of one another in» = 0, so that u,)(7) = %l- 


and this has the advantage that such a pair need only be tabulated fo 


positive values of x. For use in problems of wave propagation the mgh! 
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choice of the pair is that near infinity their oscillations should be 90° out 
of phase with one another. As will appear, these conditions imply that 
the amplitudes of either curve at its two ends are as different as possible, 
one large and the other small, and that near zero they have the character 
of the exponential function. 
For negative a the curves wu, and u, are always concave to the x-axis, 
and so resemble roughly cosine and sine functions. Nevertheless, to adopt 
the even and odd functions is again not convenient, since their asymptotic 
values again tend towards the same form. Instead it is again best to adopt 
the principle that the two asymptotic solutions near infinity should be 


those with phases 90° apart. 


2. By routine methods the even and odd solutions may be found in 
absolutely convergent series. They are normalized respectively to values 


l and 2, and are 


i | . (2 , 2s 3. 
¢ sn 50) a 
Uy - > « : ~ (\ 2a i7 ayn. (2.1) 
cK £4 n 
0 
eve i 
f i a 2 ee 
Us, \ D(gn 4 2a) re in[4)n (2.2) 
A ) 
NV Zé Ty — n! 
1 
n odd 
where n=VT4+hia), T%—=rg+ha). (2.3) 


In spite of their complex forms both these solutions are real, as is 
evident from the fact that the imaginary part of (2.1) must be a solution, 
but, unlike the real part, could only start at the term in 2? and there is 
only room for one even solution. 

Expanded for the first few terms, the solutions give 


4 76 


| ! 
i ae | ee alliage 9 
he l+a 7 (a + a1) anne (2.4) 
y3 og ; gl 
u, Y a— (a- 3) — i (a? 3a) — be cece (2.5) 
; ) 7% 5! - B: 
3. The expressions [, [D, of (2.3) play a great part in the solution. 
Write r, Ge”, lr, Gi, e%7s, (3.1) 
Also let tanB ena. (3.2) 
Then the following relations follow from the properties of gamma 
pro] 


functions 


G.G 22 cos B e774 (3.3) 


Q oan 9 
V1 Y3 p bor. (5.4) 
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Asymptotic expansions are also required. For large positive a they are: 


° 7A ? ] 5 61 
InG, = - — t1n fa +4 In 27 4+ —. 4. —__ 4. _ (3.5) 
4 . . 32a? " 256a4 1536a® 
7a l 5 61 
In G, = ——+})n}a+ 3h 227 —-—_. —- —__— — ze (3.6 
. ‘ 32a2, 256a* = =1536a 
7 7 l 7 3] 
w1+-=%3— taln ta — fa+t + —__. +} meme 4 3.7 
8 8 - ss 48a ' 5760a? ° 80640a° 


> 


For large negative a, write —a for a in (3.5) and (3.6); in (3.7) change the 
sign of 47 and replace the first term by 4aln(—4a), leaving the rest 
unaltered. 


Another relation that will be used is: 


P($+2a) = ,/(27)cos B e~74/22iaetn +79), (3.8) 


4. The expressions (2.1) and (2.2) combine very obviously into the 
solutions: 


iy = Ty y+ V2e-H AT uy = eS P(dn+ 4-4 fia)(v re! 8)"/n! (4. 
0 
Ug = T up— Vv 2e-74T wu, = ef "4# SF (—)"P(4n+ 3-4 fta)(V 2xe-7/4)" nn! (4.2) 
0 
These have the relation u,(7) = u,(—a) that is wanted, but they are 


complex. They, however, give the easiest approach for the asymptotic 
series, which can be best made by Barnes’s method. The appropriate 
integral is see 
[ D(s)P'(b—2s)(v 2ave id)2s ds, (4.3) 
ie 

taken so as to leave all the poles of the first factor on the left and of the 
second on the right. For convergence |A| < 37. By extending the integral 
round an infinite positive semicircle, over which it vanishes, it is easy to 
choose b and X so as to get a series of the form of (4.2). The result is thus: 


ix 
( ix*/4 ig ; . ; 

Us ; D(s)0(4+ 1a —28)(W 2xe-*7/4)28-1—-ta ds, (4.4 
Tri 


For the asymptotic series the Barnes contour is transferred to the left 
over a number of the poles of I'(s). The result, using (3.8), is: 


Us = 2 cos Be-87atelinit7s+7/B)y (4.5 
where 
 f ($+7a)(3+7a) . (4+7a)(8+ia)($+7a)(5-+- 1a) | 
Ua ™ = a 1 Ae 7 : < ba 


(a/v2yirial "11 (2x?) 2! (2ia?)? “y 
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5. This method cannot be applied so as to yield an asymptotic series for 
because A in (4.3) would have to be 37, which is not admissible. But 
we have as a second solution of (1.1) the conjugate complex of (4.5), and 
ean be expressed as a linear combination of these two solutions. To do 
* 


this. wu and uw, (which are real) are expressed in terms of u, and uj, and 


is given by (4.1). The result is 
Us LU , cot B + ux cosec B ¢ UY tYst7 4). (5. 1) 


))\ 


using (4.5) and (3.2). 


Us 2vrre74/46 cos B eX11t+7s-37/8)y + etml8y* (5.2) 
6. From these solutions real ones must now be constructed. They should 
have the following characters: (1) One should be the reflection of the other 
: = 0; this implies that if one is uy = Au,+ A*u¥, the other must be 
Au,+A*uz: (2) the oscillations near infinity should be 90° out of 
phase; (3) since the amplitude of uw, towards infinity is much larger than 
hat of us, they are fittingly chosen so that the product of these amplitudes 
s of order unity. These three conditions define most of the factors in A, 
uit any other real factor of order unity may be introduced. The normaliza- 
on to be used in the tables will be that the two solutions should conform 


to the equation 
du, duyy 
Uy u 


] 6.1 
‘ da dx ; ( ani 


| this will be adopted here as being best, in spite of a slight loss of 
formal elegance, since normalizations depending on squares or products 


sually, as here, introduce unwanted square roots into the formulae. 


Take A —— 5B 4/ (sec B)e™a@l4e—UM i+ ys)—ba 
iN7 . 
[he asymptotic series then are 
Uy J —t¢ “cos $f, (set 3),U, e i+Ys3)—t7 Th a 1 iar}) (6.2) 
Ur 2-te-74l28e¢ $8, (cos B){u, ea tys)+iai + yke-Uknitys+iw) (6.3) 
{gain substituting back into (4.1). (4.2). 
TF 9-in e744 (sec BIG, Uy V2G5 U4} (6.4) 
U1] 2-477—Fe7 4/4 (sec B)iG, Ug— V2G, u}. (6.5) 


hese formulae suftice to give the solutions for any value of 2, provided a 


‘not too large. They are valid whether a is positive or negative, but it 
should be noticed that whereas, when a is large positive. cos and cos 48 


ire near 1, so that (6.2) has amplitude very much greater than (6.3), when 
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a is large negative cosf tends to zero like e77, and the two asymptoti 
expressions have amplitudes of the same order. 
In the neighbourhood of the origin, the leading terms are found fro 


neh “pli ae : 
(3.5), (3.6) to be uM, 2-Ha\-i(1+a,/|a}) 16. 
Ur} 2-4 \a|-*(1—a,/ |a}) (6.7 


for large a whether positive or negative. 


7. When a is large (6.2—-5) are of little use. Take the case of a positiy: 
The convergent series (6.4), (6.5) will require a very large number of term: 
even for small values of x, and the asymptotic series (4.6) will not be usefi 
until x is considerably greater than a, whereas the curve starts oscillating 
shortly after x 2va. There is clearly need of series in which 2 and a ay 
both large, and these are provided by the well-known solution which star‘ 


i | v42?-a) dr 


e? (da?—a)!, 
[ have not succeeded in formally converting the general solution int 
one of this type, but have obtained the result by solving the differentia 
equation in this form, and then matching the result with the appropriat: 


series in (6.2), ete. This method has, of course, the weakness that there is 


no obvious remainder theorem, as there is for (4.6). 
An interesting point in the solution is that the work is made very mucl 
easier by taking the whole solution as an exponential rather than as in (4. 


as a power series. (This form is also much the best for the expression o! 


the Bessel function J, (2) when both x and n are large.) 


Let (x? — 4a) X (7.1 


: ; at+X . 
and let d?6=4| X dx = hrX—aln ; : ic 
‘ 2va 


. 
Qa 


then the solution is of the form 

’ 5 i h il m im ) a 

Inu, id—A ln X + J I eee eects Spence © i. 
4 - Bo 3 x 6 ! X 9 ! Pt 12 x 15 X 18 

where g, h, /,... are to be determined as functions of x. 


On substituting in (1.1) a sequence of equations is given, from whicl 


3g, me — 6zh, A* Ox... 


ax AL Ax 


v2 


can in turn be evaluated as polynomials in x and a. Each of these cal 


now be integrated, and in doing so the complementary functions may 
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‘Ymptotic | omitted, since they are proportional to X3, X°...., etc., and so only yield 
in over-all arbitrary constant. The results are: 
und f 
] _ ax 
q (7 4) 
( a\4s 2 
67 h (3x?+ 2a) (7.5) 
“we. 3 _,, 2, 
l | ya ax' a*z> 4+ — az3+ 19a4x (7.6) 
a?\5760 320 320 12 
siti 153 ’ . — 
m r++ 186ax*-+ 80a? (7.7) 
ot ter) bat 
be usef 1/ 3] 3] 403 4433 
scillatir n pid ax azz — a®z7° + 
: a®\ 80640 2688 2688 4032 
ind aa 
4433 ,. 13643 ,, 13049 ,., . . " 
ch star 4 asx? aa? — — a’x3— 2524a'a (7.8) 
896 SO 6 , 
6381 og, ee a! " 
p v8 29862aa4-+ 62292a227? + a? ), (7.9) 
f10N Int 4 a 
tere! T 71 
[he solution Is then 
rOpriat 
es ; iio ( ; .-% hom ) i 
tnere uu, ~ X exp/i re) g + —— pe | | = —- os Lf | (7.10) 
| X3 ao Fe" 1 ee * Fe ee i} 
ry n There is, of course, a second solution conjugate to this. 
s in (4 
SSIOI 8. The terms in g, h,... decrease in the asymptotic manner provided 
2,a, and though there is no evident remainder theorem, as there is 
" for (4.6), it may be presumed they give a good approximation. If now x 
is taken much greater than a, the series should become proportional to 
: 1.6), though the constant of proportionality might be a function of a. 
‘ ‘ Va es 
Neglecting in (7.2) d becomes 
i 
= a a 
jx?—alnx+-=mne—-, (8.1) 
and of the contributions from (7.4—9) only the first terms of g, 1, n survive, 
2: a a l 7 5) | 
whit U,—> 2 expt Ina | fe —_____ | 
Z 2° 48a © 5760a* 


80640a°| 


) 


Jial2y> j iapi{ix® +hyit+ys)} (8.2) 


by the use of (3. 


In the same approximation 
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so that as far as concerns the leading terms 


~)-—} 4 +9) 
Up, >= Z *U,erNvitys, (8.4 


[t is a remarkable fact that the asymptotic series in a on the two side 
of (8.4) should coincide so exactly, but the algebraic intricacies of th; 
Bernoulli numbers, and the complication of the solutions of the equations 
(7.4—9) are so great that a general proof of their identity for further tern 
would be difficult. If the series for uw, is expanded beyond the leading 
term in inverse powers of x, it may be verified directly that there is th 
same correspondence with wu, for the later terms, but this is unnecessary 
since both expansions satisfy the differential equation. 

For large positive a, 8 becomes nearly zero, and cosf, cos 38 becom 
unity to an approximation depending on e-7¢ which is negligible. The 
substituting for (8.4) in (6.2) and (6.3) we have 
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9. It is easy to use, the work of § 7 to give series for small x. It is onh 
necessary to write , 2 “1 
ieee ) \/(4a—2?) = iX (9.1 
and to replace id by 
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and there is a second solution, say u/., in which the sign of x is changed 
this changes the signs of %, g, /, n and leaves the rest unaltered. 

Some combination of these two solutions must be fitted so as to mate 
with (6.4) near the origin. It will suffice to compare the first two terms 
the constant term and that in x, since the satisfying of the differential 
equation will then ensure the equality of the rest. This time it will be the 
last terms in (7.4—9) that will survive, instead of the first. 

Then 
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which can be reduced to 
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by the use of (3.5), (3.6). Hence comparing with (6.4) and taking cos 8 l 


we have o uy, (9.6) 


and there is no need to take any term in w’.. Once again it is remarkable 
how exactly the two series match one another. 


10. In the case of negative a, as has been pointed out, the solutions of § 6 
stand unchanged, but when —a becomes large the approximate expressions 
are altered. The important change is that 8 approaches 90° instead of zero, 
so that cos8 = e7¢, cos$8 = 1/2 to the correct approximation. The 
changes in the asymptotic expansions of G,, ete., have been indicated in § 3. 
X is now greater than a, and ¢ must be replaced by 

x + | X dx lo X—a In x x : (10.1) 
P 24 ( a) 


0 
All the expressions (7.4-9) stand unchanged, and w, is the same as (7.10) 
with x replacing ¢. Allowing for the changed values of y, and yz the 


equation 


uy, = 2-hu, ebintys (10.2) 


still holds. Introducing the changed values of cos and cos $8 the results 
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When x is small, but a still large negative, a calculation may be made 
with these expressions similar to that of § 9, and it is found that the series 
correspond again exactly to (6.4) and (6.5). Moreover the expression (10.3) 
is also applicable for negative 2, since on changing the sign of x it passes 
over into (10.4), so that it appears to be valid for all values of x. 

It would seem that the series (10.3) may merit study by an expert 
analyst, for an asymptotic solution that is valid for both large and small, 
positive and negative values of the argument, and so also presumably for 
intermediate, is certainly remarkable. Indeed it does not fall into the 
category of asymptotic series as ordinarily defined and may well be 


convergent. The discussion would, however, be difficult because there 
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seems little likelihood of finding a general form for the terms g, h.... of thy 
series. That it is asymptotic with regard to a must be remembered, sing 
terms in e7¢ have been neglected.+ 


Conclusion 

The two standard solutions u,; and u,; adopted for equation (1.1) hay, 
the properties of being mirror images in x 0. For positive a, u; resemble: 
an exponential function near the origin; at +2va it turns and become 
oscillatory, decreasing towards +00. The amplitude is large at the positiy: 
end and small at the negative, and the phases of the oscillations are 9 
apart. For negative a the curve is oscillatory throughout, and has ampli 
tudes of equal order at the two ends. The solutions are normalized accord- 
ing to (6.1). 

The general solutions are given in (6.2—5), in which the symbols ar 
defined in §§ 2, 3, and in (4.6). 

These solutions are of little value for large a whether positive or negativ 
For positive a asymptotic series for large x are given in (8.5), (8.6), wher 
the symbols are defined in § 7. A series for small x is given in (9.3) using 
symbols defined in §§ 9 and 7. 

For large negative a the solutions are given in (10.3), (10.4) using symbols 
defined in §§ 10 and'7. These solutions are valid for both large and small 
values of x. 


+ The numerical convergence has been tested for a 10, and (10.3) is valid 
to seven places of decimals for any x. For a 10, with critical point at about 
6, there is a gap in which neither (8.5) nor (9.3) can be used ; this is from about 41 


8 for three-place accuracy, and from 3 to 9 for four-place accuracy. 
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FURTHER NOTES ON THE SOLUTION OF ALGEBRAIC 
LINEAR SIMULTANEOUS EQUATIONS 


By A. N. BLACK (Engineering Department, The University, Oxford) 


Received 28 September 1948] 


SUMMARY 
In a recent paper several methods of solving simultaneous equations and inverting 
atrices have been described. This note describes an alternative method believed 


be more convenient than any given there. 


1. Introduction 

Fox et al. (1) have described several methods for solving simultaneous 
equations, but have omitted the abbreviated Doolittle method (2, 3), 
which is shorter and more convenient than any which they give. It is 
slightly less accurate than the Choleski method, but, as will be shown 
later, this is a trifling matter. The layout given here is believed to be 


original. 


2. Method of solution 

The method used is elimination, deducing a triangular matrix from the 
given square matrix, and solving by back-substitution. But instead of 
ising the largest term as pivot, as in pivotal elimination, the terms on the 
main diagonal are used as pivots in order. This sacrifices the accuracy 
gained by using the largest term as pivot, but gains accuracy by greatly 
reducing the number of rounding-offs. 

Since the convenience of the method lies largely in the details of layout, 
the method must be described minutely. Three matrices are written down. 
A is the given matrix, of which the terms below the diagonal need not be 
written. For solution of equations the column b of absolute terms must 
be written after A. For matrix inversion no such column is needed. 
There will also be the usual sum column for checking. C is the deduced 
triangular matrix, with b and sum columns added, as for A. D is a 
triangular matrix derived from C by the rule 

d.. jit, tee >, 
/,, is not required. D has no columns for b or sum, and it must be written 
on a separate sheet, which can be folded along the lines dividing the 
columns 


‘ow | of C is the same as row 1 of A. The formation of the remaining 


; Tows proceeds cyclically as follows. Assume that r—1 rows of C and D 
: : 
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have been formed. Fold sheet D so that the rth column appears on the 
right-hand margin. Set it beside the rth and subsequent columns of ¢ 
forming successively 
r-1 
Crp = Gg + > Cyd,, for e>r. 
t=1 
Check the row sum. Form row r of D and repeat the cycle. 
When C is complete the equations can be solved by back-substitution 
It is convenient to write x as formed on a horizontal strip of paper, 
which should have —1 written in the b column. Check that the scalar 
product of x with the column sums equals the sum of the b column in the 
A matrix. In the case of matrix inversion the required b in the C matrix 
are known, for b,, = 1, and b,, = 0, for s > r. The terms 5,, for s <r 
may be non-zero, but are not required in the back-substitution. Check 
that the scalar product of each row of the inverse matrix with the row 
sums of A is equal to 1. 


3. Convenience 

The advantage of this layout is that all quantities whose product is 
required appear adjacent. The position where the answer is to be written 
can be marked by an arrow on the moving sheet. Instead of folding 
sheet D, each column in turn can be cut off as it is wanted. The single 
column is more easily slid into exact register, and the disadvantages of 
cutting up the sheet are small, as each stage is completely checked. 

The method can be readily extended to asymmetric equations, and in 
either case products are always written on the computing sheet opposite 
an arrow on the moving sheet, and dividends on the moving sheet opposite 
the divisor on the computing sheet. The sign is always that indicated bj 
the calculating machine. There is therefore little chance of error. 

An attempt has been made to compare the number of machine opera- 
tions in this method with that of Choleski, for the case of six equations 
Such comparisons cannot be definite, for they depend in part on the 
machine used and the user’s technique; for instance on whether successive 
divisions by the same divisor are done by direct division or multiplication 
by the reciprocal, and in the latter case whether the reciprocal is formed 
by break-down or build-up division. 

The table below shows the comparison between the two methods, both 
for solution of equations and for inversion. Setting means setting a numbe! 
on keys or levers, for summing, multiplication, or division. Summing 
the operation of turning a number, already set, once into the product 
register. Square rooting means the complete operation of obtaining 4 
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square root; this has not been broken down into its elements as the 
process depends on the number of figures needed. 


ne set of equations Vatrix inversion 
’ holeskt Doolittle Choleskt Doolittle 
Set . , 216 218 329 279 
Sur ; 100 III 108 87 
[ 8 86 1809 156 
28 21 37 36 
S ° 6 


[It will be seen that the balance of advantage lies slightly but definitely 
with the Doolittle method for solution of one set of equations. This 


advantage is greater if the complete inversion is required. 


4, Accuracy 

The solution by this method of the set of equations treated by Fox et al. 
1) is given below. The errors in the fifth decimal of x (cf. p. 162) are 28, 
15, 51, 27, 29, 15, so that the accuracy is superior to their elimination 
process, equal to that of Morris, and inferior to the other methods, Choleski 
in particular. However, these errors correspond to an error in the caleu- 
lated coefticients of 2 to 3 units of the last decimal carried in the working. 
They are therefore fully covered by the retention of one place of decimals 
is guarding figure, which would be required in any case. Any labour 
expended in attempting greater accuracy is wasted for the reasons stated 


in the paper under discussion 








= 
x Xs Y pa 
27647200 + 18469100 |+ 1236790 + 244615400 
18469100 28026900 |+ 4844800 (+ 271177701 
26397400 + 16793600 |+ 12495006 + 20403320 
16793600 + 26324600 47 3040 + 210136000 
20157800 11492100 1064700 + 131604200 
19406500 375 310C 124295900 
157200 120512500 48508 20 
Cc 
27647200 158469100 1236790¢ + 244615400 
8322418 + 10129420 7140313 |+ 34133101 
2301425 4060659 |+ 953036 + 13944342 
793359 31727% $62700 |- 314745 
765461 |+ 1316771 583101 |+ 2665330 
63536 267067 |— 203528 
197¢ 51198613 |— 34202100 
7 29712339 - 30103019 
7594513 172196807 
73718549 
x 
7 799258435 420339650 100000000 
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Although it does not arise with this set of equations, it should he 
emphasized that accuracy cannot be obtained unless all the terms on the 
principal diagonal and the largest term'in each b column are brought to 
the same order of size by suitable changes in the roots of the equations, 
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THE WAVE RESISTANCE OF A CYLINDER 
STARTED FROM REST 

By T. H. HAVELOCK (King’s College, Newcastle-on-T yne) 

[Received 13 August 1948] 


SUMMARY 


A method of obtaining expressions for wave resistance in accelerated motion is 
given, but the particular problem examined is the motion due to a circular cylinder 
submerged at a given depth below the free surface, the cylinder being suddenly 
started from rest and made to move with uniform velocity. The surface elevation 
it any time is discussed, and expressions obtained for the finite wave resistance 
at any time after the start. Numerical calculations have been made for three 
lifferent speeds, and curves are given showing how the resistance rises initially 


and oscillates about the steady value for each speed. 


1. Introduction 

CALCULATIONS of wave resistance have hitherto been made only for a body 
moving with constant velocity, the problem being treated directly as one 
of a steady state when referred to axes moving with the body. The case 
of non-uniform motion is of interest in itself, and also has possible applica- 
tions. For instance, in measuring the resistance of ship models, the question 
arises how long it is before the effect of the starting conditions becomes 
inappreciable. As a matter of fact, measured resistance curves always 
show oscillations about the steady value for a given speed, but these are 
no doubt mainly due to the natural period of the measuring apparatus; 
however, it would be of interest to have some examination of the approach 
to the steady resistance after the initial stage of accelerated motion. 

Expressions for wave resistance in accelerated motion have been given 
by Sretensky (1), who obtained them by transforming the hydrodynamical 
equations to a form suitable for axes moving with acceleration, but the 
formulae are too complicated for numerical calculations in general; 
Sretensky has, it is understood, made some calculations more recently for 
a particular law of acceleration, but the results are not available. 

In some early work (2), instead of assuming the steady state as estab- 
lished, I used an alternative method for uniform motion. This may be 
described as finding the disturbance due to an infinitesimal step in the 
motion of the be dy and then integrating. It was pointed out at the time 
that the method could be applied to motion with variable velocity. It is 
shown now that this method leads directly to expressions equivalent to 
those obtained otherwise by Sretensky. However, the present note deals 
mainly with one particular problem, namely, a circular cylinder submerged 
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in water at a given depth, suddenly started from rest with a given velocity 
and maintained at that speed. It has been found possible to make numerical 
calculations in this case, and the results illustrate various points of interest. 


2. Circular cylinder 

Take the origin O at the centre of the circular section, of radius a, at 
a depth f below the free surface, with Ox horizontal and Oy vertically 
upwards. If the cylinder is given a small horizontal displacement ¢ 67, from 
rest to rest, the velocity potential of the subsequent fluid motion is given 
by m 
b = 2ca*g! Sr | e-*?/-Y sin(xx)sin(g!te!)x? de. (1) 

0 

This is equation (12) of the paper already quoted (2), obtained there by 
a Fourier integral method; it can also be derived in the manner given later 
by Lamb (3) for the three-dimensional case. 

The velocity potential for continuous motion with variable velocity can 
be found by a direct integration. We consider first the simple case when 
the cylinder is suddenly started at time ¢ = 0 and made to move with 
uniform velocity c. We obtain, noting that the origin is at the centre of 
the moving cylinder, 

4 cate — cate : 

x+y x*-+-(27—y)* 
f $ 
+2ca*gt | dr | e-*°S-Y” sin{x(a+ct—cr)}sin{g'«!(t—7)}«! dx. (2) 


0 0 


Deriving the surface elevation » from the relation 





CY Ch , 
I — (y =f). (3) 
ct cy 
we obtain 
t L 
n 2ca® | dr | e-*S sin{x(a+ct—cr)cos{g'«!(t—7)}« dx. (4) 
0 0 
Hence we have 
2a*f , f cosxx 
1] “> Ts nga? e-"S dx 
a-++-F* J K—Ko 
0 
x 
fe a ee en ee ee 
~ ca" pene = = . seh Nas "Sdx, (8) 
| Ke+g'x? Ke—gi! 


0 


where x, = g/c?, and the principal values of the integrals are to be taken. 
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WAVE RESISTANCE OF A CYLINDER 


The first two terms in (5) give 


af ’ . ae 
' —_ — LK, A°E- “0! SIN Ky X — 
c 1-7 
>» [| Kysin«f—Kcoskf 
K“-- Ko 
0 
2a*f ) 
"1 a. a7 Ky, A" SIN Ky X — 
yet f2 
, 45 [ kgSin«f—Kcoskf R 
2x, a? y——— etd (x <0). (6) 
J K“-+- Ko 
0 


The expressions for », represent a steady state relative to the moving 
cylinder and symmetrical fore and aft of it. With 2+ct = € = distance 
from a fixed origin at the starting-point, the last two terms in (5) represent 
the surface elevation at any time due to an initial displacement and velocity 
which is the negative of that given by (6); this must be the case in the 
present problem and it can be directly verified. With a change of variable, 
and with w2 = x,y, the last two terms in (5) are given by the real part of 

 pilgut—git . piléu?+gtu) 5 
” 2a? ure" du—4a? | - ure du. (7) 
u—u J Utuy 


0 


The limiting value of », as t becomes infinite is derived from the principal 


value of the first integral in (7); taking the real part, we find that 


Ko] 
7 


lo > 2K, a" snk,x as t>-+0o. (8) 

Turning to (6), we see that ultimately (8) cancels out the regular waves 
in advance of the cylinder and doubles the amplitude of those in the rear. 
Without examining the surface elevation in detail, we may specify more 
closely the part which at any time consists of a regular train of waves 
accompanying the moving cylinder. It is clear from the form of the integrals 
in (7) that the only contribution to such a train comes from the first 


integral, or from 


: p eiléu? g?tu) ae 
2a2u2 et du, (9) 
p u Ug 
and is due to the pole at uv = uw». Regarding u as a complex variable, the 
path is along the real axis indented at u = uy. There is a saddle-point at 
w= gt 2. First suppose € > 0. The path of integration may be rotated 


round the saddle-point to the line of steepest descent, namely, the line 
u = gt 2€+re*'”, the contribution of the circular ares required to complete 
the closed contour being zero in the limit. We have also to take account 
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of the indentation at uw) according as uy > or < g't/2€, that is, according 
as € > or < $ct. In this manner, it is found that as far as the regula 
waves are concerned, (9) gives 
2rK,are-“ sinkyx (E > $ct), 
-2rK, ae sinkyx (€ < det). (10 


Similarly if € < 0 the line of steepest descent is the line vu = g*t/2€+-reli= 
and the corresponding contribution is — 27K, a?e-"°! sin Ky x. 

Summing up this outline of an analysis, the surface elevation at any 
time is made up of three parts: (i) the local symmetrical disturbance 
travelling with the cylinder given by the first and third terms in (6): 
(ii) a regular train of waves 47x, a7e-*/ sin x, 2 behind the cylinder extend- 
ing from x = Otox Sct; (iii) the part given by the remaining integrals 
representing a disturbance which spreads out in both directions and 
diminishes in magnitude as time goes on. 

The second part agrees with the general description using the idea of 
group velocity. The third part has not been examined in detail, but an 
asymptotic expansion suitable for large values of € and t may be found 
from the transformed integrals indicated in the previous discussion. For 
large positive values of € and (g't/2é+)—u, &?, the first term in such an 


expansion is vag! 
mighat? 


2€4(¢— ct) 


For € = ct, that is, at a point over the centre of the moving cylinder 


eI cos(4ta—gt? 4¢). (1 


this reduces to 

at(“0)*« ‘xo! cos }(7—Ky ct), (12 

. ct 

a result which can be obtained directly from the integrals in (7) by using 
the method of stationary phase. After a sufficient time, (12) gives approxi- 
mately the departure of the motion over the cylinder from the quasi- 
steady state consisting of the local symmetrical disturbance and the regular 
train of waves to the rear. If A, (= 27/x,) is the wave-length in the regular 
train, the wave-length of the disturbance near the cylinder is 4A,, the wave- 
length for which c is the group velocity. The usual direct solution for motion 
with uniform velocity leads to the surface elevation (6) with regular waves 
in advance as well as to the rear. The so-called practical solution is then 
obtained by superposing a free infinite wave train cancelling out those in 
advance and doubling the amplitude to the rear. Another well-known 
method of obtaining this practical solution directly is to use the frictiona! 
coefficient introduced by Rayleigh. In the present analysis we have not 
used this frictional coefficient, the values of the integrals being interpreted 
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as principal values wherever necessary. The chief point of the discussion 
is that there is no ambiguity when the motion starts from rest. The motion 
which is gradually established as time goes on is the practical solution for 
the steady state, with regular waves only to the rear of the cylinder; this 
result is in fact associated with the group velocity being less than the wave 


velocity. 


3. The wave resistance 

The wave resistance may now be obtained to the same approximation 
as for steady motion. From (2), withz = x+iyandV = (g/x)*, the complex 
potential function is 
ca ca" . ‘ ‘ . > — 

re ll cag dr f ix(e—V\\t—r)__ p—ix(c+V Ml Te —iKz 2h eb dx. 
0 0 (13) 
In this case the resistance is given by the usual form of Blasius’s theorem, 


namely, the real part of 


,_: [ (dw\? 

Loi dz, (14) 
7” aw dz 

taken round a small contour enclosing the origin. Expanding (13) in the 

neighbourhood of the origin, this gives 


; 


R 2mpc*giat ( dr ( f ix(c—V \t—7) e—in(e+V yl Micke 2«h dic 
0 0 
Sorg pK5 a -~ prer2kol” dy, (15) 
| v—l v+1 | 


U0 


where the real part is to be taken. The value to which R tends as t becomes 
infinite arises from the first part of the integral; and we find that 


R—- 4r?qpn2ate-2F as t>to, 16 
YP*0 


which is the resistance for uniform velocity. 
With « K, ct, B 2x, f, the form 


‘ * (sin av(v—1 sin av(v-+1)) - 
R SorgpKs, a4 = z : x )\ 15 


e-Be* dy, 17 
v—] v+i1 | a m4) 


0 
is suitable for numerical computation by direct quadrature if « is not too 
large and 3 not too small. 


For other values of the parameters it is possible to find more suitable 
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forms by transforming the integrals in (15). We divide the integrals into 
two parts. Put 


r (y>— ] ; ae a ee ‘ 
I, | eater __ ciae| piar®—Br? dys (18) 
Jot 
0 
x 1 
and kL | cfav*—Br* sin av dv = k | e-tentd—0") dig, (19) 
ry 0 


with p (B+ta)/(B?+ 7), and k Lap. 

The second form for L in (19) can be evaluated by quadrature, or from 
convergent or asymptotic series which can readily be found. It can be 
shown that 


i, = G+ —zy A+ it say ® + 
X a 
| 1—ki—k? +k) + k4+-(1 — 3ki—ok2)* +. 3 kL. (20) 
| ¥ v7| 
The remaining integral from (15) is 
. iar(v—1) iar(v+1) 
L= ft i | e-Br? dy, (21) 
* J | el v-+l1 | 
0 


Considering v as a complex variable, the path for the first term in (21) 
is along the real axis indented at the point v = 1. We take it now along 
the real axis from v = 0 to v = } and then along the line v = $+,re?’, 
with 7 positive; in the limit, the contribution of the circular are required 
to complete the contour is zero. The second term is taken along the real 
axis from v = 0 to v = —}3, and then along the line v —3+rei'7, with 
r positive. In this way we obtain 

4 


> 


i mie-B-+2 | eiar(v—1)—Be? qa (v—1)+ 


0 
‘ 
. f fe Biret'7 +43)? p—B(reti7_4y° 
xX)e 


e-or dr. (22 


atin OS + 


\ reti7_} 


| 
reine 
0 
The integrals in (22) are not too difficult for computation for fairly large 
values of « and moderate values of 8. The resistance R is given by the real 
part of —4zgpx5a*i(Z,+ J). 
When a is sufficiently large, the divergence of the resistance from the 
steady value becomes of the order a~?; this first term in an asymptoti 
expansion comes from the first term in (18) and the second term in (22). 
Applying the method of stationary phase to these integrals leads to the 
result 


R — 4rg px? at{re-2*0f — 2 (a/x, ct)'e-!*0f sin }(7—K,yct)}. (29 
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WAVE RESISTANCE OF A CYLINDER 331 
The second term in (23) corresponds to the approximate perturbation 
in the surface elevation given in (12); this can be checked by direct calcula- 
tion of the extra resistance due to (12). 
4. Calculations of the wave resistance 
The steady resistance (16) is a maximum at the speed for which 
3—2«,f— 2. Calculations have been made for the resistance in the 
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present problem for three different speeds, given by Bf = 2, 4, 6. The 
expression (17) was found suitable for values of « up to about 15, taking 
in interval of 0-05 for v and not attempting any high degree of accuracy. 
The forms (20) and (22) were used for higher values; these and (17) being 
calculated for some values in common for checking purposes. Finally, (23) 
was found sufficient for very large values of a. About fifteen spots were 
calculated at each speed, and the results are sho.’n in curves of resistance. 

The abscissae are values of €/A,, where € is the distance travelled from 
the start and A 
ordinates are proportional to the wave-resistance R, the broken straight 


in 
» is the wave-length for the corresponding speed. The 
lines showing the steady resistance at each speed; the wave-lengths for 
the three speeds are 2zf, xf, and 43af. 

The following seem to be the chief points of interest. Considering that 
the steady resistance is associated with the idea of group velocity and the 
establishment of a train of waves to the rear, the resistance rises to the 
neighbourhood of the steady value more rapidly than might have been 








332 T. H. HAVELOCK 


expected, especially at the higher speeds. Another point is that the maxi- 
mum resistance attained is larger relatively at the lower speeds; thus at 
the highest speed shown the maximum is about 1-11 times the steady value 
at that speed, while for the lowest speed the maximum is about twice the 
steady value. After about ten wave-lengths the resistance is given approxi- 
mately by (23) and the amplitude of the oscillations diminishes slowly 
with the distance travelled; but even at an earlier stage the period of the 
oscillations is roughly 4A). 


5. Three-dimensional motion 

Turning to the general problem, we take the origin O in the free surface 
with Oz vertically upwards. Suppose a point source of strength m is 
suddenly created at the point (0,0, —f) and maintained for a short time 
dr. To satisfy the condition at the free surface for the initial motion, we 
take 


$= =", (2 


with r= = x?7+y?+(z+f)*?; 72 = a? +y?+(2z-f). 

The initial surface velocity found from (24) acting for a time 87 gives 
a surface elevation which can be put in the form 
ba or | dé | e-*f cos(kar)K dk, (25 

7 0 

with wo = xcosé+ysin 8. 

The velocity potential of the fluid motion at any subsequent time ¢ due 
to this initial displacement without velocity is 


d g* dr | dé | e-*I+*= Cos(«ar)sin(gitic! i? dk. (26 

0 
Consider now a source moving parallel to Oz at constant depth f, the 
strength m being a function of the time. Let 2 be measured from a moving 
origin vertically over the source, é from a fixed origin at the starting-point 
and let s, be the €-coordinate of the source at any time ¢. Then we obtain 


from (26), 
t 7 


mit mt) gf i . ‘ eee ae 
re) ) Lg m(r)dr | dé e-*I+* cos(Kar’ )sinfgi«!(t—r) x? ds 
4. rs wT. d ‘ 
0 T 0 (27 
with ow’ = (€— $,)cos 6+ y sin é. 


We may generalize this result for a solid body moving through the liquid 


If the solid moves through an infinite liquid with unit velocity, we maj 
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take the fluid motion to be that due to a certain distribution of sources 
and sinks over its surface and of amount o per unit area at each point. We 
assume this distribution in the present problem in order to obtain the wave 
motion to the first approximation. Thus in (27) we replace m by ac(t), 
where c is the velocity at time t; if (hk, —f) is a point on the surface of 
the body we also put «—h for x, y—k for y, and 
Coa (€—h—s_,)cos 0+ (y—k)sin 8, 
and the required velocity potential is obtained by integrating over the 
surface of the body. 
We shall not carry the general problem further meantime, but consider 
the case of a slender ship form. Here the usual approximation is to take 
(ey/eh)/27, where the surface of the form is given as an equation 
for y in terms of h and /; further, the source distribution is taken to be in 
the longitudinal section of the form by the plane y = 0. We obtain, in 


- 
this case, 


eff) [fit l \ oy ; 

d ~ dhdf 

Ln | | (- ro) oh ; 
: [ { “anap { e(r) dr ao [ KI+*2 eos(Ka’ )sinfgix?(t—r) x! dx 
—s5 ; dhdf | c(t) da | at é cos(Ka@’ )sinig?«?(t—T) }«? dk, 
is ; 7 6 (28) 


with ow’ (€—h—s.)cosé+ysin@. This result is equivalent to that 
obtained by Sretensky by a different method. 
The pressure at any point is given by p é¢/ét, neglecting the square of the 


fluid velocity; and the resistance is found from 
- i ; vin im —— 
R 21 | pr’.o, 4 an'as’, (29) 
2 ch 
taken over the longitudinal vertical section. Hence, from (28) and (29), 
we find 


R | away 





[ { (aay + (fN(e aye + (+f | dha 4 
Ja , ; dh 
ad t - 
gp cy ; ; red y . aos 
! “3 | | zich df | | ah dhdf c(t) dr dé » 
| e-+Neos(Kar’ Joos{gta'(t—r)}x de, (30) 


0 


with w’ (h’ h + S; 8_)cos @. 


‘ 
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The coefficient of ¢ is an effective mass for this particular problem, 
taking account of the free surface and assuming no wave formation and 
noting that the square of the fluid velocity has been neglected. 

As a special case, suppose the model to be started from rest with a 
velocity c which is then maintained constant. The finite resistance at any 
time after the start is given by the second term of (30), with ¢ a constant 


and 
a’ = fh’—h+c(t—r)}cos 0. 


The result can be reduced to the form 


t 7 


R gpe | dr dé (12+ J?)cosf{xe(t—7)cos Oeosfg'x!(t—r)\« dk, 
"se ow & 31 
VW ith I 1 iJ : | | Cy P —~Kf+ikh cos 9 dhdf. (32 
JJ eh ’ 


Integrating with respect to 7, this gives 


p= 20 faa f esa 


0 0 


sin{(xec cos0—g'x!)t! | sinf(«e cos 8+-g'x!)t} 


KC cos 6 —g?k? Kc cos 6+-g' x! 





It can be verified readily that the limiting value to which this tends as 
t becomes infinite is 


le 


R = 2X09? | (72+ J2)sec36 dd, (34 
8 
where J,+-7J, is given by (32) with « replaced by x,sec*é. 
This result (34) is the known expression for the steady resistance at 
constant speed. 


REFERENCES 
1. L. SkRETENSKY, Joukovsky Cent. Inst. Rep. 319 (Moscow, 1937). 
. T. H. Have tock, Proc. Roy. Soc. A, 93 (1917), 520. 
3. H. Lams, Proc. Roy. Soc. A, 111 (1926), 14. 


no 








ON 


By I 


Th 
of th 
harm 
sym! 
whic 
to th 
moti 
studi 
depe 


is un 


4 
IN a 
cire 
posi 
can 
pro 
prir 
em} 
the: 
oTe 


a PR 


Ye 


har 


eXi 


vith 


iT any 


nstant 


nee 














ON THE ROLLING MOTION OF CYLINDERS IN THE 
SURFACE OF A FLUID?+ 
By F. URSELL (Department of Mathematics, The University, Manchester) 
Received 14 October 1948] 


SUMMARY 
The methods of a previous paper (ref. 1) are extended to permit the evaluation 
of the fluid motion when a cylinder of arbitrary symmetrical section performs simple 
harmonic oscillations of small amplitude about a mean position with its axis of 


symmetry level with the surface. Special attention is paid to the slow motions, for 


which a method of successive approximation is developed. Applications are made 
to the virtual mass in slow heaving, to the determination of the roll axis, and to 


motion in a sea-way. Finally, the waves generated by the forced rolling motion are 


studied. It is shown that their amplitude, and thus also the damping decrement, 
depends critically on the details of the section. A section is calculated whose rolling 


s undamped to the first order. Experiments are needed. 


1. Introduction 

[x a former paper (ref. 1) the motion of a fluid was considered when a long 
circular cylinder executes vertical heaving oscillations about a mean 
position in the surface. It will be shown that the method of that paper 
can be extended to discuss the oscillations of a long cylinder of any section, 
provided that the cylinder is symmetrical about a vertical plane and the 
principal axes of the section are of the same order of magnitude. The 
emphasis will be generally on skew-symmetrical (rolling) motions. In 
these the length of the waves produced by the rolling is usually much 
greater than the beam of the ship. Application of a method of successive 


approximation outlined in ref. 1 is therefore permissible. 


2. Formulation of the problem 
On the assumption that the motion is two-dimensional and simple 
harmonie of period 27/c, a velocity potential ¢ and a stream-function 


exist satisfying 2 


= = 0, (1) 


at = 0, (2) 
Ci coy” 
where the origin of coordinates is at the mean position of the centre of the 
cylinder and y increases with depth, the mean free surface being y = 0, as 
in ref. 1. The stream function is prescribed on the cylinder, or, to a first 
Read to the VIIth International Congress of Applied Mechanics. 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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approximation, on its mean position. To the same order, the condition for 


the free surface is ad 
Kéd+—=0 (y=0) 9 

oy 
where gk = o*. (4 


Introduce dimensionless coordinates £, 7 defined by the isogonal trans. 
formation 


x 
a =ccoshésinn+e > Ag,,,¢-@r+M sin(2r+ 1)n, (5 
r=] 


y =csinhécosy—c > A, 
1 


or+1 & -2r +E eos(2r +-1)y. (6 
Fe 


The constants A are chosen in such a way that the given cylinder 


2r+1 


becomes the curve € = &, and the region occupied by the fluid is mapped 


isogonally on —iz7 < yn < da, & < &€ < «©. Equation (1) becomes 
- : 2 ] 37, So ¢ | 
oo a9 
ak 7h 
—st+a5 = 9 (7 
0€* § an? 


and the surface condition becomes 
r e 2 9 né od 1 
Ke|sinh € + SY (—1)r-te-2r+ fA +3] +—=0 (7 +4, & > &). 
a > on 
(3) 
3. Skew-symmetric motions 
Consider in the first place, motions in which 
$(€,n) = —4(€, —7n), (9 
for instance a rolling motion. It is sufficient to consider the region 
0 < » <4. Laplace’s equation (7) and the surface condition are both 
satisfied by the skew-symmetric set 
-donsalés) = €-2"40E sin(2n-+1)q— 
nan ¢8in(2n+2)y _ 
4n 4n+4 


sin2n7 , 
ee adh: 


ie Ke| € 2n€é 


x“ 


a +] In | Oy 19 
Ss Ae ,, etansarsng ant ort 2)n 
. a 2n+2r+2 
- 


(n = 1,2,3,...). (10) 


? 
It will often be convenient to work with the conjugate stream-functions 
YonailE, 7) e-(2n+DE eos(2n + 1)» 


7) ae 
_ Ke| e-2ng0O8 227, o-c2n ,9¢C08(2n-+-2)n _ 
| 
4n 4n+4 


oa pictus) Oo. 9 
_ r ; (2n+2r+2¢C08(2N-+ 2r-+ 2)n 
; : et Maal Tn 
‘ wae In+2r+2 
‘ ; 


(nm = 1,2,3,,...). (11) 
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ion for [It is not to be expected that the general skew-symmetric motion can be 
expanded in terms of these sets, since at infinity the motion consists of 
regular wave trains travelling away from the origin, whereas each 
donat(é> 7) V anishes at infinity. 

t The set is completed by the addition of the potential ® and stream- 
trans- | function ‘’ describing a horizontal dipole at the origin (cf. ref. 2). In 
particular, the stream-function ’ is given by 


Y ¥(E, 7)cos ot 4 V(é 7 )sin ot 


j 


9 q 

ae 7 f kee-*= dk 

Im | 2K cie**2-‘ot 1 “sin ot —— 
k+iK 


ylinder 
: a ZC 
lappe Im | 2K cie*ke—tot 4 -S1n of - 





2 = “ (iKzy .,.. — (iKz)/l | 
‘a <a S ot) ik Cly log il 2 % Mt — I » = oe 
- _— sie heal r! boii La =o! +r ; 


\ 


r=0 r=] 

4. Rolling oscillations of small amplitude about the origin 
The process which is employed to determine the coefficients of y,,,,(€, 7), 

8 n= 1,2,3,... in the expansion of % is best illustrated by an example. 

Consider a cylinder whose equation is € = €, and suppose that it is rolling 

about the origin, the angular displacement at time ¢ being 

Q 6 4, cos ot, (13) 


where 6, is small 


The boundary condition is (ref. 3, Art. 72) 

ub s(x" y”)+ Fit), (14) 

on the cylinder, where F(t) is a function of t only; for small 6, this con- 
dition may be replaced by 

1 dé 


yp r2ty*)+ F(t) on &€ = &. (15) 
For other skew-symmetrical motions x?+-y* is replaced by some other 
even function of y, to which all the following considerations may be 
»t10nNs ° , 2 ° ° 
_ applied.) Suppose that % is multiplied by a constant so as to make the 


coefficient of ’ unity; A is a function of 0, and the period. Then, 


Aus YY" (§ cos ot ‘Y.(€, 1)sin ot- 


Porsi( Ke, €)for+1(€, 7) C08 ot - 
1 


Jor +4 (Ke, Eo)ho,.4(€, 7)sin at, (16) 


(11 r=] 


Z 
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and P41; Yori, must be determined in such a way that on € = &, the 
stream-function ¢% differs from a multiple of x?+-y? by a function of ¢ only 


that is ‘E.(Eg, n)cos ot ++ '¥,(E,, 7) sin ot — 


Por sa(Ke, Fo) oy+1(Eo, 9) C08 ot — 


| 
Mes 


| 
“Mes 


Jor+4( Ke, Eq)tbo,+1(Eq, )8in ot 


= O(t)(2*+-y")o+ F(t). (17 

There must be no source-singularities on the cylinder, i.e. ys is continuous 

in the range 0 < 7 < 3m. It is therefore assumed that the series (16 
converges uniformly in 7 on € = &}. 

Suppose that the beam of the cylinder is 2a and that the draught is } 


Put 7 = $7 and subtract the new equation from (17), thus eliminating F(t), 


[ Va(Eo. 7) —‘F.(Eo, $77) |cos ot +-['F,(Eq, 1) — Fe (& , $77) |sin ot — 


Por+1lPorsa(So> 1) —Yor4a(€g, $77) ]oos ot — 


“Ms 


Jor tl Yor+a(€o. 7) —Por+1(€os jm) |sin ot 


| 
“Ms 


- O(t)(a?+-y?—a"), 


> 


:, xv?+y*?—a?’ 
(p, cos of +-q, sin ot)( a) » say; (18 
0 


that is 


¥(y, 2) —Ve(Eos 77) 


¥ (Eo; ”) E3(€o, 377) 

v2 + y?— a? a 
= as( = a? “) v 3 Jor+1|Pors1(Eos 2) —YPorsa(Eo,37)]- (20 

0 1 
These equations must be solved for the infinite set of unknown 
Porst> Yors4: One approximate method is to fit a finite set of polynomial: 
by least squares, as was done for the case of the heaving circular cylinde 
in ref. 1. However, in this paper the solution of the set of equations wil 
be considered only under the condition that Kee is small and that a an 

6 are of the same order of magnitude. 

The expansion in power series of Ke of the functions ¥,(£, 7), Ygé.1 
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suggests that 


= 3) 
Por+1 > Py or+a( Ke)” (21) 
n=0 
x o 
% (1) F a (2) ¥ 9¢ 
Yor +] > Qn er 1\ Ke r T log Ke > mer+1( 48)", (22) 
n=0 n=1 


since the expression in square brackets on the right-hand side of (19) and 

20) involves only (Kc)° and (Ke)*. It is therefore possible to determine the 
tn,2r tn,2r 

(Ke)" and (Ke)"log Ke on the two sides of (19) and (20). It is assumed 

that the resulting power series defining p9,,1, Yo,,, converge for sufficiently 

small values of Kc. This is known to be true for the heaving motion of the 


coefficients P, o.4, Qd-.1, Ol... recursively by equating coefficients of 


circular cylinder. 


Now, if (Kc)*log Ke is negligible, 


W(é,7) = 2Ke, (23) 
mM. 2se — 
¥(E, 7) in| —~+—i1Ke(y- logike)], (24) 
V.(Eo, 1) —¥lEo 4x) = 0. (25) 
; _ 2c : | Z| 
W5(€o, 7) — ‘¥e(€o: $77) Im (=) +=Kelog|* 
Teo 7 1a | 
2c . 
Im/. ) + (Ke), (26) 
7% 
and bora (Eq: 1) —Wop aa (Eo, 477) = er +DSe cos(2r+- 1)y+O(Ke). (27) 
It follows that — Phorsa = Pharr = QPhrar = 0 (28) 
and (53.4, is given by 
2, (2 4 y2— q?\ dl 
[m| -| Qu | — ] + D> Oba e-@r+ocos(2r+1)y. (29) 
Tg) ” a 0 —_ 
rT 


The amplitude of the motion of the cylinder is determined to the first 
order by Qj}, which may be found by multiplying (29) by cosy and 
integrating from 0 to 47. But 

Z¢ —< - a 
Im (—  ° Aone" *V§ cos(2n+-1)n (€ > &)), (30) 
from the theory of harmonic functions, where only a, is required in the 
calculation of QM 


0.1" 
Multiply both sides of (30) by e€ and let € tend to infinity. For large €, 


a ~ 4ce§ sin », (31) 
y ~ 4ceé cos 7. (32) 
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In the limit the equation (30) reduces to 
—(4/7) cos yn = a, cos n, 





or a, = —4/7. (33) 
On carrying out the multiplication and integration, 
$7 
4 4 { /(x?+y?—a? 
——e-f = Qi) — | ea ) cos 7 dx. (34) 
T 7 | r—a* Jo 


0 


On the cylinder, the stream-function is given approximately by 


9 9 9 
a?-y?—a? , a i . 
Ay = QO sin a(S) —4A6,osin ot(x?+y?—a*),, (35) 

_ 0 
where A is a normalizing constant, and at positive infinity 
, , . gAB_ , . 
Ais = 2Kce- ®¥ cos(Ka—ot) = J: —e- Kv eos(Ka—ot), (36) 
Co 


where | B| is the amplitude at infinity. Hence 
K%b(0*—a") _retn [ sos oa , 
= — oa — = Kcek§, (a*— x —y”), cos 7 dy. (37) 


B 





. 
0 


As an example consider the elliptic cylinder 


x = ccosh & sin », y = csinh &, cos n 
a = ccosh &, b = csinh &. (38) 
Here (a®—a*—y?), = (a®—b?)cos?n (39) 
and the amplitude |B} is 
$77 
Kce&f), ( |a®—b?\cos*n dn = 3K?0)(a+b)?|a—b]. (40) 


0 
[t can be shown that this formula is valid whether a is greater or less 
than b. The calculation shows that the error is O( K3a‘), which is negligible 
if Ka is small. When Q§}} has been found, Q0},,, is obtained by multiplying 
(29) by cos(27+1)y and integrating from 0 to 47. 


5. The boundary condition é¢/éy = 0 (y = 0) 
There is another way of approaching the problem, which is suggested 
by the boundary condition 


it =0 (y= 0). (41 
cy 
As K tends to zero, the boundary condition tends formally to 
op =0 (y 0). (42) 


oy 
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It is not unreasonable to expect a close relationship between the first 
approximation to the wave problem and the complete solution %, of the 
problem with boundary condition (42). (Cf. ref. 10.) One way of solving 
the latter is to expand the stream-function in terms of the set 


é 


(22 +E cos(2n+- 1)», (43) 
to which is added the dipole function 


Tr im(=). (44) 


TT 


Then the determination of the coefficients in this expansion is clearly 
the same as the determination of Q$2,,,.. Expansion of Im(2c/7z) leads to 
the following 


THEOREM 1. Suppose that in the skew-symmetric wave-problem x is pre- 
scribed (except for an additive constant) on the cylinder £ = &, and that at 
the mean free surface y 0 ; ed 
K d -— = 0. 

oy 

T'o obtain the wave amplitude at infinity when Ka is small, solve the simpler 
problem, in which yb is prescribed as before on & = &, but the free surface is 
replaced by a rigid plane and 


S—0 (y=0). 
oy 
If the solution to this latter problem is written in the form 


wr > Baas’ (2r+DE cos(2Qr+ 1)7 sin ot, 
r=0 


where B,, 1é (2r+1fo (5). (45) 
r2 


Then for the original problem the waves at infinity are described by 
us bor B, Kee-*®¥ cos(K \x|—ot). (46) 
The error in xs is small of order K*a®B, when B, + 0. 
A similar investigation for symmetrical motions, using a source function 
and a set of polynomials y,,(€, 7) (ef. ref. 1), leads to 
THEOREM 2. Jf the solution of the problem modified as in Theorem 1 is 


written in the form 


2r 


x 
br , { 077 > sf 


ert sin 2rq|sin ot, (47) 
1 


where Cy 1s adjuste d to make 
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then for the original problem the waves at infinity are described by 
wb = nC,e-*¥ sin( K |x|—ot), (48) 


the error in being of order K*a*C, log Ka, when Cy + 0. 


6. The pressure on the cylinder 


The pressure at any point in the fluid is connected with the velocity 
potential by the relation 


od 
Pp = 9py—p— + Fit), 


where F(t) is a function of time to be determined from the condition that 
at the free surface the pressure is zero. In the skew-symmetric motion ¢ 
is considered as the superposition of the non-orthogonal polynomials, 
defined by equation (10), and the dipole potential taken in the form 
Sa 
® = Real part of | 2K cie'* tot 4 sin at sa a" : 


0 





(49) 


The non-orthogonal set tends to zero as |z| tends to infinity, (49) reduces 
rt tho; 1 set tends t tends to infinity, (49) red 
to the potential describing a regular wave train for which the pressure is 
given by ob , 
P = 9pY— P= = 0 on the surface. (50) 

C 


Hence without loss of generality F(t) may be taken as zero. A simple 
calculation gives 
THEOREM 3. If in the modified skew-symmetrical motion 


[epp/ey = 9, y= 0| 


Up > Bops. 6-2 cos(2r+ 1)y sin ot, 
r=0 
2r+é F 
where By, ete = Of}, (51 
7s 

, - ; : ed 
then in the original motion the hydrodynamic pressure Ps near the 

Cc 


cylinder can be obtained approximately from 


x 


dp = — > Boye? sin(2r-+1)nsinot. (52) 
. 


The error in the pressure is of order Kcp @¢,/ét. A rather more complicated 
calculation shows that if the modified symmetrical motion is given by 


2 
re 


~ . + 9 ] 69 
br [° ont > Cy, e-** sin 2rn| cosot, where C,,e-*"% = O (" (53) 
r=1 
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¢ 


, oD . ‘ 
then the hydrodynamic pressure “Pa, +8 obtained from 
Cc 


r 


oo 
bn = [—Colytlog4Ke+é)+ ¥ OC, 
r=1 


e-*rf cos 2rn| cos ot+-mC, sin ot. (54) 
7. Applications 

The theory which has been given in the previous sections may be 
applied to a variety of problems concerning the motion of cylinders, and 
under suitable conditions to the discussion of the motion of ships of long 
parallel middle body. If numerical computation is to be avoided, the 
parameter Ka must be so small that the first approximation is valid. 
This condition is usually fulfilled when a long ship is given a slow 
forced rolling motion, but not for a heaving motion. Also, the error in 
the calculation of a symmetrical motion is of order Kalog Ka, and is 
larger than for a skew-symmetrical one. For these reasons most of the 
calculations will be made for cases of rolling motions. The section of the 
cylinder will first be taken as elliptical, so that direct comparison with 
experiments on ships will not be possible. The extension to other sections 
involves rather tiresome calculations which, though simple in principle, 
would obscure the ideas underlying the work. However, the calculation of 
the amplitude of waves generated by the rolling of a cylinder of nearly 
rectangular section will be given in detail. 

For convenience the modified motion in which the free surface is supposed 

cal od , 
to be replaced by a rigid plate | -2O,9 0) will be described as the 

oy 

R-motion set up by the motion of the cylinder. 

A) The roll axis of a cylinde r of elliptic section 

When a ship executes a slow rolling motion in still water, there is a 
point of least movement, called the tranquil point, about which the rolling 
motion may be supposed to take place. The position of the roll axis in 
the corresponding two-dimensional motion will now be investigated. 

Consider a ship the immersed part ot which is a cylinder of semi-elliptic 
section with its centre in the mean free surface, and consider the skew- 
symmetrical forced rolling motion about any axis parallel to the axis of 
the cylinder and lying in the plane of symmetry. For slow rolling the 
first approximation is adequate; also, the damping is so small that during 
a roll period there is little difference between free and forced rolling. To 
maintain a motion at constant amplitude a small amount of energy must 
be supplied per unit time, most conveniently by a couple. In general 
there will be an oscillatory reaction on the roll axis, but when this axis 
coincides with the axis of free rolling, the reaction will be small. For small 
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angles of roll, the rolling motion about a fixed axis can be considered ag 
made up of a rolling motion about the centre, together with an oscillatory 
horizontal displacement of the centre in phase with the rolling motion, 
As the cylinder rolls, the fluid exerts a varying pressure on the cylinder, 
given by ad 
P = 9PY— P=» 
ot 
the first term representing the hydrostatic, the second the hydrodynamic 
pressure. The potential is adjusted in accordance with Theorem 3. By 
the principle of Archimedes the total hydrostatic force is vertical and 
depends only on the displacement, which is constant (to the first order 
in 9,) during the rolling motion and equal to the weight of the ship. 
The horizontal forces acting on the cylinder are: 


(55 


(1) the resultant force due to the hydrodynamic pressure; 

(2) the reaction of the roll axis. 

The resultant of these two forces is equal to the product of the ship's 
mass and the horizontal acceleration of the centre of gravity. When the 
roll axis is the axis of the free rolling motion, the second force will be 
neglected. An equation is thus obtained to determine the roll axis. 

Let d be the depth of the centre of gravity and / that of the roll axis 
below the water-line. ‘It is assumed that the horizontal axis is the major 
axis. Suppose that the inclination of the ship at time ¢ is 

6 = @ sinot, (56) 
then the horizontal displacement of the centre of the ellipse is 
—l§ = —l@)sin ot. 
It follows, as in ref. (3), that on the ellipse € = & 
ys = —4}6,0c*(1+cos 27)cos ot —10, oc sinh &, cos 7 cos at 
se r—1 
=— a) ac? cos of pa ei arb3) cos(2r+-1)y— 
ro 
—l@,o0csinh €,cosycosat. (57) 
The solution for the corresponding R-motion is clearly 


~ 





4 “3 (—1)2 | 
— —6, ac? cos ot > a e-2r+Xé-L eog(2r-+-1)n— 
—16, oc sinh &, e—€-cos 7 cos at. (58) 


By Theorem 3, the corresponding potential function is 


4 ; = (—1)-1 
_= bs) 2 Ss t 
dr - 9 GC" COS a i (Qr— 1)(2r-+1)(2r-+-3) 


+10, oc sinh £,e-¢-9sin y cosot —(59) 


e~2r+IX$—L0)gin (27x +1) q+ 





near the cylinder. 
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The horizontal force per unit length opposing the motion is 


$7r $7 
— , ’ Oh Oy 
— | pay = 2 | Ady & =o 
J ov ot on 
— $7 0 
pco sinh &,| 349 oc? +47, ocl sinh £,|sin ot. (60) 


The horizontal displacement of the centre of gravity is 
(d—1)@, sin ot. 
Now the horizontal force per unit length is equal to the product of the 
mass M per unit length and the horizontal acceleration of the centre 
of gravity: 
Mo*(d—1)0, sin ot = peo sinh €, sin of[ 36) oc? +- 476, ocl sinh & | (61) 
i.e. (with ccosh 5 = a, csinh &, = 6), 
M(d—l) pb[ $c?+- 4b], 
_ Md—% pbc? 


so {= —___., 62 
M +-42pb? (62) 
But the cylinder is in equilibrium in its mean position, so that 
M = }xpab, 
_45(q2— ph? 
1 — ad sm(a b ) (63) 
a+b 


To this approximation the position of the roll axis is independent of the 
frequency. It follows that (63) also gives the position of the roll axis for 
the case of free slow rolling, which can be expressed as the superposition 
of a number of undamped forced oscillations with periods lying in a 
narrow band near the mean rolling period. Equation (62) was first derived 
by R. Brard (ref. 10), who compared it with experiments and found 
satisfactory agreement. 

When the semi-axes a and 6 are nearly equal, 


1 = 4d. (64) 


This formula is in agreement with the experimental results mentioned 
by Sir William White (ref. 8, p. 169) for rolling without bilge keels. It 
appears that the position of the roll-axis depends critically on the differ- 
ence (a—b) of the semi-axes. When higher terms are included in the 
calculation, it is no longer permissible to assume that the two skew- 
symmetrical motions are in phase. Under these conditions, i.e. when the 
diameter of the ellipse is not small compared with the wave-length, there 
is no roll-axis; the instantaneous centre of rotation is not confined to a 
small region 
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(B) The added mass in a slow heaving motion 
Suppose that an elliptic cylinder is given a forced heaving motion 
y = lsinot. 


Then on the cylinder é = &, 





ys = —loc cosh €,sin 7 cos ot 

locoosh §, [294-2 5 {— 
= —loccosh &, | - — > ——~—— sin 2rn | coset. 65 
tolc0t5 2 nae ”| ‘ (66) 

It follows that in the corresponding R-motion 

~f2.,2— (-1y7 ' 
wp = —loc cosh £,|—»+ - >; e-*rE-f0sin 2rn |cosat. (66) 
a" Ly r(4=—1) 


By Theorem 3, the corresponding potential function ¢ is 


lac co hé Sf loo 1Ke €) | - : (—iy* e 2r(f—Eq) ) | t 
— bloc COSD ¢9| ——(y- 108 AC) --— >) oes 2—2nS—S0 cos 27'n | cos ot — 
7 , r(4r*—1) 








77 


—2loc cosh &,sin at. (67 


The vertical hydrodynamic force opposing the motion of the cylinder is 


7 $7 
| p—dn = —2p | a “dy (€ = &), (68 
J on ot On 
-—to7 0 
9 ¥ 
= 2pce cosh £, lo*c cosh £9 —(y+log $Kc-+ &5) | cos n dy sin ot— 


0 


47 $77 


ioe) 
-})r-1 , 
on | cos 2rn cos 7 dy sin ot +2 | Cos 7 dy coset 
, (472 
1 
0 


0 


oof 2/. 1. Keef\ . 2. 1< l 
2 pla?o? [= (y-Hoe —)sin ot——sin of Ware t 7008 at (69 
1 
es) 2/1 ; af 
= 2pla*o”| ——| log —_—_ + 0-23 } sin ot +-2 cos ot}, (70 
a\ °K(a+b) 


where the numerical constant 0-23 = 1-5—log 2—y approximately, and y 
is Euler’s constant. 

Here the first term is 180 degrees out of phase with the acceleration and 
so represents the effect of added mass. The second term is in quadrature 
with the acceleration; the work done by it is accounted for by wave 
damping. 
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Let M’ be the mass per unit length of a semicircular cylinder of radius 


a, i.e. M’ = }rpa?. 


Then the inertia force per unit length represented by the first term can 
be written in the form 








Mt’ |log—_* 4-0-3 |<¥, (71) 
17 —K(a+b) | dt? 
whence the added mass per unit length is to the first order 
l 
=’ log y ! 0-23 e 72 
| K(a+b)* (72) 


It may be noted that the ratio of the inertia force (71) to the hydrostatic 
restoring force per unit length 2palg is of order Kalog Ka, which is 
assumed to be small. 

(C) Motion of a ship in a sea-way 

Suppose that a ship of the type described in (A) is placed in a progressive 
train of waves advancing from positive infinity 
gA ee ( K A oe 
ob = ——e-KY sin( Ka-+ ot) J" — (x cos ot —y sin ot)+- G(t) (73) 

o Co 
near the origin to the first order); and suppose that as a result the 
coordinates of the centre of the cylinder at time ¢ are 
« = Lcos(ot+-«), y = lsin(ot+-f), 


and that the inclination of the minor axis at time t is 


G 6, cos(at 4-8). 
Then, on the ellipse & é,, 
% = —loc cosh €, sin y cos(ot-+-8)— Loc sinh €, cos 7 sin(ot + «)+- 
145 oc? sin(ot-+ 5)| 1-—-cos 2). (74) 


The stream-function satisfying this condition is 


le a os 
e-Kv sin( Ka-+ot) 
qA a P y 9 1)’ 1 . : 
— Ke cosh &, cos ot | — y I 2 | ‘ y e— 2-0 gin 2rn — 
o T 7 “x r(4r2 -1) 
| 2 2— | ee 
— loc cosh €, COS(at >) | => ae are -£)sin 2ry a 
| 7 7 a 1) 
L i : 
gA ,, . . (t—£,) . ' \ 
Ke sinh €, e~-*% cos n sin ot — Loc sinh €, e~¢-*% cos n sin(ot+-«)+ 
4 : “ ( om rl ° ¢ \ 
+-@,0c?sin(ot-+-8) S —_ e~2r+ ES cos(2r+-1)n. 
- LZ, (2r—1)(2r+1)(2r+3) 


(75) 
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Whence the potential function is to the first order 


$= e. 


Ky cos(Ka-+-ot)— 


— a=" cos ot +loa cos(at +B)|| —=(y-+og$ +6) + 
Co 


2~— (—1y-* wee 
t= > _ e-2nt-Ldeos 2 
7 & r(4r2—1)\ _ r|— 








_2 PS sin ot +loa sin(ot + p)| ne 
o 
Kt i ; £ . . 
or : 2 ~(€-S gin 7 sin ot + Lobe-€-S0sin 7 sin(ot--a)— 
oc 
4 © Ir | 
— —4 ac? sin( (ot + 5) > a XE gin(Qr + 1 Ye. 
- a ee r+1)(2r+3) 


(76 
From this velocity potential valid near the cylinder the force and the 
couple acting on the cylinder can be found, the expression for the pressure 
being 
5 C d 
P = 9PY— P=: 
ot 


When the forces are equated to the product of mass and acceleration in 
the direction of the force, these equations, with a similar regiage'e for the 
couple, are sufficient to determine the quantities 1, LD, 0); «, B, 8. For 
example, consider the heaving motion which is sain of the skew- 
symmetrical motions. 

The force opposing the motion is (correct to the second order) 


anesthe 
1 
~~ ] + 0-23 || A sin ot +1 sin(ot+- 
= pa® o ee ai)" I sin ot +1 sin(ot+-8)]+ 


+ £ pa?o*[ A cos ot +-1 cos(ot-+-8)|+ 2pgla sin(at-+-8), 
7 


which is to be equated to 4mplabo* sin(ot+-B) (77 


since the mass of the ship is equal to }zpab by the principle of Archimedes. 
It follows that to the second order 
Bp = 0, l+A=0 (78) 
(taking account of the equation gK = o?), so that the ship moves with 
the wave. 
The reflection from the ship’s side due to heaving can be shown to be 
zero, to this order. 
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8. The problem of roll resistance in still water 

The damping of the rolling motion is known to be mainly due to three 
different causes: (1) skin friction, (2) wave-making, (3) eddy-making by 
keels. It is agreed that skin friction does not account for more than a 
small fraction of the rolling decrement (ref. 2). The waves generated in 
the rolling motion of a ship are of the order of one inch in height and 
rather difficult to measure. Nevertheless, naval architects have arrived 
at tentative conclusions about the relative importance of (2) and (3). 
For instance, G. 8S. Baker has developed an approximate theory (ref. 4) 
from which he deduces that wave-making should account in some cases 
for more than half the energy lost by the ship. According to Baker, (2) 
and (3) should make comparable contributions to the damping for any 
section of rectangular shape with rounded corners. The contribution from 
eddies is increased when the ship is under way. 

The mechanism of wave-damping is independent of the viscosity of the 
fluid. A reasonable estimate of its magnitude may probably be obtained 
by assuming that the fluid is frictionless. If the ship has a long parallel 
middle body of length at least as great as the wave-length set up by the 
rolling motion, the two-dimensional theory worked out in this paper may 
be applied. Moreover, the parameter Ka is small; in most cases its value 
is about 0-1 (ref. 2). It is therefore sufficient to consider the first approxi- 
mation. The calculation is unchanged when a uniform flow parallel to 
the axis of the cylinder is superposed. 

[t will now be shown that hydrodynamical theory does not support the 
empirical formulae hitherto put forward. In the next section the waves 
generated by a cylinder rolling about a point in the mean surface will be 
calculated. (The tranquil point is usually near the mean surface.) It will 
be shown that there is at least one nearly rectangular section for which 
the wave amplitude is zero, so that the damping may be overestimated 
by Baker’s formulae. It is hoped to arrange experiments to check the 
theory. 


9. Determination of the waves made by a conventional section 
The formula (37) for the wave amplitude requires that 2 and y should 
be of the form (5, 6), with € = &,. 
Transformations are known expressing a rectangle with rounded corners 
in this form (ref. 5), but for practical purposes it is necessary to use finite 
expressions 


2N 

Lv(C, Eq) ccosh €)sin y+-¢ > Ao.44 ye-"tesin(2r+ 1)n, (79) 
1 
2N 

Yn(C, Eq csinh €, COS n—C > Aeris went Dg cos(2r-+- ] )n; (80) 
. , 
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for which x?,+-y- can be found without a prohibitive amount of computa. 

NTYN puta 
tion. Equations (79, 80) will represent a roughly rectangular curve jf 
contact with its tangents at » = 0 and » = 3m is of the highest order 





possible: 9 \2s 
(=) YN = v, n = Q (s = a (81 
C 7) 
(=) ty=0, n=4r (s=1,2,...,N). (82 
07) 
Che equations for A,,,, y are 
sinh €, = > A oy43,n C+ D50( 274-1) (¢ = 1,2,...,N), (83 
T= 
2N P 
cosh €, Dd Aapasry en erthso(—1)r-1(2r+1)* (s = 1,2,...,N). (84 
r=1 


By addition and subtraction, 


N 
efo — > Ay,_, e-4r-Mho(4r—1)% (¢ = 1,2....,N), (85 
r=1 
N 
—}e-te = > A gp Cr TSO 4 4 Ly (¢ = I,3....,N). (86 
r=1 


Each of these systems is of the form 


N 
> B,ybf=1 (8 =0,1,...,N—1) (87 
r=1 


the solution of which is 





= (b,—1 . 
aii [| ( 5} o 


The equations (83, 84) thus have the solutions 


(4s—1)?—1 (89 
ifs 1)2—(4r. nh | 


‘ 1)?—] 
A —(4r+D§0( 47 t-])2 etl aan o (90) 
—— ail la i] | Fearn (474 ” 


As N tends to infinity each coefficient mae to a limit corresponding t 


A yoy (4r—2So( 47 — 1)? pete | 


a 





the limiting rectangle. Small values of N give satisfactory sections, somé 
of which are shown in Fig. 1. The parameter €, determines the shape 
the parameter c the scale. It is necessary to verify in each case that ther 
is a (1,1) correspondence between points on the boundary and thi 
coordinates € or 7. 

In the subsequent calculation N will be taken as unity. Then thi 
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parametric equations of the cylinder are 

a = $(a+6)sin 7 +3(a—b)sin n+ 3(a+5)sin 3n—Z(a—b)sin 5yn, (91) 

) 


25(a —b)cos n—j(a+b)cos 3n+4(a—b)cos 5n, (92 


y = $(a+b)cos n—® 


also cege 8(a+b). 


The amplitude at infinity has been shown to be, by (37), 
K*6, ces | (a?—ax?—y?)cos 7 dy). 
0 
In this case 


) 


g2--y2—a? 


67 


128 


(a*—b?)(1-+-cos 2n)-+ 
-i:5| 25(a—b)?+ 81(a-+-b)*|(1—cos 4n) +-33;(a®—b?)(1+-cos6m), (93) 


whence 





7 
bs (x?+-y?—a*)cos 7 dy a (a? b*) + — [25(a - -b)?+-81(a+b)?| 
7. : L052 2707 
. 9.9 
= — ~~" (a—1-26b)(a+1-05b) (94) 


and the amplitude at infinity is, to the first order, 
0:63.K76,(a+b)(a+1-05b) |a— 1-266). (95) 
When a and 6 are equal the amplitude is 
2729 @3 ‘ 
1 K*6,a%. (96) 
This may be compared with 2K26,a* for a thin plate (ref. 9). But when 
a = 1-26b, the first order amplitude vanishes. It may be deduced that 
for this ratio of a to b, the amplitude, and therefore also the wave-damping, 
is very small. The corresponding section is shown in Fig. 1. 
For reference the amplitude is given for the same class of cylinders 
turning about an axis at height H above the water surface. It is 


0-63K26,(a+-b)| (1-26b—a)(1-05b+-a)+-0-015H (a+ 26b)]. (97) 


10. Discussion 

That wave-making might account for the damping of the rolling motion 
of ships was first suggested by W. Froude (ref. 7). Assuming that the 
angular displacement satisfied a differential equation of the second degree 
with constant coefficients, he was able to give the damping decrement 
in terms of the characteristics of the ship and the wave height at infinity 
(see also ref. 2). A considerable advance was made by Havelock (ref. 2), 
who not only corrected an error in Froude’s calculation but also attempted 


to find the wave height at infinity in terms of the characteristics of the 
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ship. For this purpose he considered the motion as two-dimensional and 
evaluated the waves made by the rolling about a vertical mean position 
of a completely submerged thin plate whose width was small compared 
with the length of surface waves it generated. Comparison with Baker’s 
experiments showed that the calculated amplitude was too small. This 
was to be expected, since there is a flow between the surface and a sub- 
merged cylinder, which tends to reduce the amplitude, and which cannot 
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Fic. 1. Some typical sections. 
Parametric equations: 


a : : 
I. 2= ‘. (9 sin 7+sin 3y), 


a rf ~~ . - 
y= z (9cos y—cos 3y). (Minimum radius of curvature 0-25a.) 
f 


49a / 3 : 
IL. 2 = Fy, (15sin 7+2sin 39— 7 sin 79), 
t= 7g (15 sin sin 39 — 75 sin 79 
49a “ F 3 ~ = . . > 
Y = S40 { 15 cos y—2 cos 3n+ Fo °° in). (Minimum radius of curvature 0-16a. 
9 , 25 , ] ] — 
Lil. «= 16 (e+ b)sin 7-4 ia (* —b)sin n+ 16 (4+ b)sin 39— Za (4 -b)sin 5n, 
9 25 : | Z 
16 (a- b )jcos ams 48 (a b)cos n 7] G2 Lh jcos 3 — 438 '7 —b)ecos oO; 
where a = 1-266. When this section is forced to roll about the origin the wave amplitude 


of the resulting fluid motion vanishes to the first order. 


take place when the cylinder is in the surface. This effect was calculated 
by the present writer, who gave an explicit expression for the waves 
generated by a thin plate oscillating about the vertical and not completely 
submerged. Comparison with experiment showed excellent agreement 


(ref. 9). 
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In those calculations it was assumed that it would be sufficiently 
accurate to replace a cylinder by a thin plate of width equal to the 
draught of the ship. When the finite section of the cylinder is taken into 
account, formula (95) is obtained. Baker’s experimental results, as quoted 
by Havelock, belong to a ship in which a = 1-316 in the centre section 
(ref. 4), so that the good agreement with the earlier formula cannot be 
explained on any two-dimensional theory. Dr. Baker has informed the 
present writer that the model ship on which he made experiments did not 
possess a long parallel middle body, so that two-dimensional theory could 
not be expected to apply. No detailed comparison of theory and experi- 
ment is possible until experiments are made on suitable models. 
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THE YAWED DELTA WING AT INCIDENCE 
AT SUPERSONIC SPEEDS 


By GWENDOLEN M. ROPER (Imperial College, London) 
[Received 27 April 1948; revised 8 February 1949] 
SUMMARY 

Three cases of the yawed delta wing at incidence have been considered: (1) the 
case when both leading edges lie inside the Mach cone of the vertex; (2) the case 
when both leading edges lie outside the Mach cone of the vertex; (3) the case when 
one leading edge lies inside, and one leading edge lies outside, the Mach cone of the 
vertex. In each case the wing has two leading edges, and the trailing edge lies 
outside the Mach cones of the base vertices. The formulae for the disturbanc 
velocities in the regions outside and inside the Mach cone of the vertex have been 
obtained. Hence the pressure distribution on the isosceles delta wing has been 
deduced for each case, and the total lift calculated. Some numerical results ar 
exhibited graphically. 


I. Introduction 

1. In a previous paper (1) we have considered the flat delta wing at 
incidence at supersonic speeds, when the leading edges lie outside the 
Mach cone of the vertex (hereinafter termed simply ‘the Mach cone’), the 
wing being symmetric to the direction of flow. 

The present paper deals with three cases of the yawed delta wing, with 
the leading edges outside or inside the Mach cone. 

The analysis for the case when both leading edges lie within the Mach 
cone is reduced, by a suitable transformation, to that for the symmetri 
delta wing with the leading edges within the Mach cone given by H. J 
Stewart (2). 

For the cases when one or both leading edges lie outside the Mach cone 
the methods of the complex variable are used. to obtain the induced 
velocity components inside the Mach cone. The conditions outside and 
on the Mach cone are quoted from known solutions. 

For each case, the pressure distribution on the wing and the total lift 
are found, 


II. Formulation of the problems 

2. In each case the x-axis is taken through the vertex of the wing 
parallel to the undisturbed stream, the axes of y and z to complete 4 
right-handed orthogonal system, the y-axis being in the plane of the wing. 
The speed of the undisturbed stream is U, its Mach number M, and the 
Mach angle » = cosec-1M. The angles between the leading edges of the 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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wing and the x-axis are w,, —w,, and the (small) angle of incidence is a; 
in applying the boundary conditions, z is put equal to zero. 

The velocity components of the disturbance in the x, y, z directions 
respectively are denoted by u, v, w. The thickness of the wing is ignored. 
The equations used are those used for the symmetric case, (2), (1), that 


is. the Prandtl-Glauert equation of the linear perturbation theory 


,2P @P ep 


9 ‘ * 9 ‘ 9 
Ox* Oy* ——s G2* 





; (2.1) 
where P is any one of the velocity components u, v, w, and B denotes 


(M2—1)!, and the irrotational relations 


ow ov ou ow ov ou 


= - — Sa — (2 2) 
' —, —. 2.2 
oy Oz oz Ox Gx oy 
The condition to be satisfied on the wing, that is on 
Z 0, £ Q, -xtanw, <y < rtana,, 
is that tan « w/(U+-u), or to our order of approximation, 
w= —Ua. (2.3) 


3. In the previous paper on the flat delta wing (1) it was deduced from 
(2.2) that the equation to be satisfied by the induced 


equations (2.1), 
velocity components, inside the Mach cone, is 


eP @P 
a = 0, 3.1 
Co- oe? ( 
where Q(o/2_1. »2\8 
sech o mt ed Btanw, | 
x (3.2) 
6 = tan-1(z/y), 


and that the equation to be satisfied by the induced velocity components, 


outside the Mach cone, is 


2p 2p 
A (3.3) 
ox” ob? 
where 
sec x 3(Yy* 2°) . Btan w, | (3 4) 
6 tan-!(z/y). 


w, # are the normal spherical polar angular coordinates, o is real within 
the Mach cone, and is zero on the Mach cone, x is real outside the Mach 


cone, and is zero on the Mach cone. 
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Solutions of the equation (3.1), for the region inside the Mach cone. 
will now be obtained for the yawed delta wing for the cases (w, > w,): 


(a) 0< w, <p, 0 < wy < p, that is, both leading edges within the 
Mach cone. 

(b) wy > p, w. > p, that is, both leading edges outside the Mach cone. 

(c) w, > p, 0O< wy < p, w,—w, < 7—2p, that is, one leading edge 
outside the Mach cone, and one leading edge inside the Mach cone, 
and the trailing edge outside the Mach cones of the base vertices. 


Solutions of equation (3.3) can be found and the constant velocities 
outside and on the Mach cone deduced mathematically; or it can be 
argued, on physical grounds, that the conditions correspond to those for 
the infinite yawed wing. Since these conditions are well known, the 
results will be quoted. 


Ill. The yawed delta wing, when the leading edges of the wing lie 
within the Mach cone of the vertex 
4. For the region inside the Mach cone the velocity components u, v, v 
satisfy the equation 


eP ep 
——eefoa——- x @, (3.1) 
Co~ 0G? 

Therefore, the complex functions U = u+ia, V = v+i6, W = w+ii, are 


taken as analytic functions of the complex variable # = 6+-ic. It has 
been shown (1) that 


dV dw 
- cot d (4.1) 
dd dé 
dU dW 
and —— = —— O0800 F ——. (4.2 
dd B dé 
Inside the Mach cone, 0 < o < «0, on the Mach cone, o = 0, —7 < 6 <7, 
and on the leading edges of the wing, 
secho = secho, = Btanw, (@ = 0), 
secho = secho, = Btanw, (60 = 7). 
Therefore, on the wing, 
6 0, %<0< 0, 
or @=4+7, o<ac<@. 
On the #-plane, the Mach cone is represented by the 6-axis, from @ = —7 
to @ = +7, and the upper surface of the wing by the two semi-infinite 


lines 0 = 0 (0, Ko <0), 0=7 (0, <0 <@). 
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The transformation 


{= £+in = d+A/(f+cos#), 


> > i 


(4.3) 


vhere —d, —f, A are positive constants, and f < 1, transforms the semi- 











nfinite region bounded by the lines 6=0, 0=aw (¢ > 0); o=0 
)<6<7) in the #-plane, into the region » > 0 in the ¢-plane. The 
. 00 F Aco 
E fo, 
Bo; 
- 11 0 —_— 
C D @ 
Fic. 1. #-plane. The points A, B, C, D, E, F in Fig. 1 correspond 


to points A, B, C, D, E, F in Fig. 2. 


| boundary is mapped on the whole €-axis. The values of d, f, A are chosen 


so that, in the ¢-plane, the points corresponding to points (7,0), (7, ¢), 
0,c,), (0,0) in the #-plane, are the points —1, —k, k, 1 on the €-axis. 
It is easy to show that these values of d, f, A are given by: 








f d sinh 3(¢,—«a,)cosech 3(a,-++a,), 
, ee a (4.4) 
A 1—f* = sinho, sinh o, cosech? $(¢,+-0,), j 
AN 
- 00 -| -k 0 +k +i +00 
+ ---- rn ——e > 
» € : B OC 3 
Fic. 2. {-plane. The points A, B, C, D, E, F in Fig. 2 correspond 
to points A, B, C, D, E, F in Fig. 1. 
and therefore k cosh }(¢,—o,)sech $(0,+0,). (4.5) 


The transformation 
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that is € = —ns’, modulus k, maps the region » > 0, in the ¢-plane, into 
the rectangle having its corners at points (’ = +K, +K-+iK’ on the 
¢’-plane, K, 1K’ being the quarter-periods of the Jacobian elliptic func- 
tions, modulus k. The points ¢ = —1, —k, k, 1, © in the ¢-plane corre- 
spond to the points K, K+7K’, —K+ik’, —K, 0 respectively in the 
¢’-plane. 

H. J. Stewart (2) has shown that the conditions for W are satisfied by 





aad = iDcd?(¢’), D being a constant, (4.7) 

C 

where —Ua = Re (Mar, (4.8) 
0 

and therefore D= - a (4.9) 


where E’ is the complete elliptic integral of the second kind, with modulus 
k’ = (1—k?)*. Therefore 





dW dW dt’ a k?Ua (f2?— ma (4.10) 
at xa 806lh ee... k?) 
Hence, using equation: (4.2), 
dU _ tne me nD (4.11) 
a ~ BEA Py\(C—e) | 
and J Ue { = f_| (4.12) 


«BE UP \(C—B (C—B) 
There is no integration constant, since, on the Mach cone, o = 0, u = 0. 
On the upper surface of the wing, tanw = y/x, secho = |Btanw!, @=0 
or 7, n = 0, and 


_atfeosho+1  f+Ptanw 


t - 4.13) 
6 = bo +eosho+f  1+4/8tanw ( 
Therefore, on the wing, on the upper surface, 
oe Ux { K®—f&o \ (4.14) 


~ BEF? \ (ey 
When w, = @, 0. = 0,,f = 0,k = secho, = Btanw,, and on the wing, 
= fBtanw, and therefore 
== = U x tan? a — (4.15) 
E’(tan’w, —tan*w)! 
which is Stewart’s result for the symmetric wing. 
From equation (4.1), 


dV ik?Ua { fl—1 


a ~~ Fa PNe—eyy 
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, Ux {( C—fk* | 
>, INto n | hence Vy _- = (4.16) 
ant 7 > 9 > 
= tie Bf? (eB 
func. | There is no integration constant, since v = 0 on the Mach cone. 
corre When w, = w,, on the upper surface of the wing, 
n the , —Uatanw 
E’(tan’w, —tan®w)*’ 
ed by | which is also the value obtained by integrating Stewart’s expression for 
(4.7 dV in the symmetric case. 
5. The lift of the yawed delta wing, when both leading edges lie within the 
(4.8 Mach cone of the vertex 
If the disturbance velocities wu, v, w are small compared with the un- 
as disturbed stream velocity U’, and their first powers only are retained, the 
(4.9 ; : : ; 9 
pressure on a surface element of the wing may be calculated by Bernoulli’s 
| equation . 
dulus | dp 1a oe 
L oq constant, 
i a 
(4.10) and we obtain p—p, = Ap pUu, p being the density of the un- 
disturbed stream. 
Hence, on the upper surface of the wing 
BE'(1—f?)* | (k2—&)4J 
4 19) Since the pressure variation on the under surface of the wing has the 
opposite sign, the pressure difference is 2|Ap|. 
0. To find the total lift, the pressure difference is integrated over the surface 
7—() of the isosceles delta wing. On the wing, R? = 2?+-y?, tanw = y/x and 
c is the maximum chord of the wing, measured perpendicular to the 
413 trailing edge. The angle between the chord c and the z-axis is 
A (w,—Wo)/2, (5.2) 
and the total lift is 
I W ¢ sec(w—A) 
4.14 2pU? * k?—fé F 
;) a a = Jo dw RdR 
BE’(1—f?)# (k2—€?) . 
wing, a= 
y k 9 
pUncI—f yh f_U=feddée (a. 
7 9 9 9? Vert 
4.15 E ‘, (k?—&)*(H, €)—,)? 
where H, BfcosA—sindA, H, BcosA—f sina. 
Hence, the total lift is 
Uac?(1—f?)tk?2 (H,—H. . 
pU “ac*(1—f*) (415 if) (5.4) 





E’ (H3—H? k?)i’ 
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whilst the lift coefficient based on area is 


CO. = 27raBk?(1—f?)!cosA - 
L ~ Fi’ tan }(w,-+o,)(H2— Hk)! ar 
ma cos A 


= [cos 2A—cos 2y cos 2y— 
E’ cos uw cos y ~ 


— 2{sin(u+y-++A)sin(u+-y—A)sin(u—y-+A)sin(u—y—a)}}}t 

(5.6) 

where the semi-angle of the wing is y = (w,+w,)/2, and yu is the Mach 
angle. 


When w, = w, the total lift is 


mpU*axc* tan*w, 27a tan w, 


Re and C, = Pa A (5.7 
which are the results obtained by H. J. Stewart for the symmetric case, 
When wa, = p, o, = 0,& = 1, f —1, so that 
1—f?=0 and H3—H?k? = 0. 
1—f2 1 
But lim — = =tan(u—A 
dim (asa) = pene) 


and therefore 
; 4xcosAf1—BtanA]} . 
C.> sony 4a cos A[tan(u—A)tan p}}, (5.8) 
Bt B+-tandA 
which gives the lift coefficient for the delta wing with one leading edge 
within the Mach cone, and one leading edge in the Mach cone. This value 
agrees with that given later in (11.7). 


IV. The yawed delta wing, with both leading edges outside the 
Mach cone of the vertex 
6. Solution for the region outside the Mach cone 
The equation to be satisfied by the velocity components u, v, w for the 
region outside the Mach cone is 
Pp _&P 0. (3.3) 
ex? o6? 


x is real outside the Mach cone, and is zero on the Mach cone, and 


sec x = |Btanw|. 


On the leading edges, 


sec y = secy, = Btanw, (@= 0), 


(6.1) 
sec y = secy, = Btanw, (0 = 7). 
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It can be shown that the values of u, v, w for the regions outside, and on, 


the Mach cone are given by: 


0 0 ex u Uy 2 vy 
t—(X2—X) 0 T u Us a —Vs, \w —Ux. 
—(¥,— X) 0 0, u —, 2 ~) (6.2) 
T 6 [7 (Xo x)| u —Us 2 Vo, 
‘=F 0 T—(¥e— ¥); u 0, y= 9, w = 0, 
where TT /(8 8} ; Une 
Uy x/(P SIN x1), Vv; —Uacot x;, | (6.3) 
Us Ua (6 sin x5), Vs —Ux«cot X2: | 





Fic. 3. AOB, HOF are the Mach planes of the leading edges OA, OH of the 
wing. From the geometry of the figure, the angles subtended at the axis OC, 
of the Mach cone of the vertex, by the Mach planes AOB, HOF, are given by 


1/B tan p tanw,cos ACB tanw,cosHCF. 


Therefore angle ACB = y,, and angle HCF = y,, since B tanw, = sec y,, and 
Btanw, = secy,. Angle ACD = y,— x and angle HCE = x,— x. 


7. Solution for the region inside the Mach cone 

Che equations used for the solution inside the Mach cone are equations 
(3.1), (4.1), (4.2). 

Inside the Mach cone, 0 < o <0, on the Mach cone, o = 0, —7 <<0< 7; 
and on the wing. @ 0 or +7. 


On the #-plane the Mach cone is represented by the 6-axis from 6 = —7 
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to @ = +7, and the upper surface of the wing inside the Mach cone, by 
the semi-infinite lines @ = 0, 0 = z (o > 0). 
Asoo sg 


-TT -(1-X,) -X; 0 X, T-X2 T 
B C D E 3 
Fic. 4, #-plane. The points A, B, C, D, E, F in Fig. 4 correspond 
to points A, B, C, D, E, F in Fig. 5. 











The transformation C= é+in —cost (7.1) 
transforms the semi-infinite region bounded by the lines 6 = 0, @=7 
(¢ > 0); o= 0 (0 < 6 < 7) in the #-plane, into the region » > 0 in the 
¢-plane. The boundary is mapped on the whole €-axis. The points —1, 
—cos x,, +cos x, +1 on the €-axis correspond to the points 0, Xp 7X» 
7 on the @-axis. 





a1 
-00 “| ~COS X} 0 cosx, | +00 
A B C DE FE 





Fic. 5. ¢-plane. The points A, B, C, D, E, F in Fig. 5 correspond to points 
A, B, C, D, E, F in Fig. 4. 


The conditions to be satisfied by w inside the Mach cone are 


w= WwW, —Ua 
on the wing, and 
0< 0 Xr» | w= w - 
o 0 =< 
™—Xo < |0 a7, | on the Mach cone. (7.2 


X1 < |0| < 7— x2, w= V. 


Therefore, for the region considered here, the conditions to be satisfied are: 


— .< one w= —Ua 
cos x, < § < +0, } on the é-axis. (7.3) 
—cos x, < COS Xo, 
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There is a discontinuity in the value of w when € = —cos xy, and when 


& COS Xo- 


Conditions (7.3) are satisfied by 


W = w+ *og( oe), (7.4) 
{+cos x, 


the principal values of the logarithms being taken. Therefore, on the 


upper surface of the wing, 


: Wt —cos #?— cos y, — 
W oe loe( x3) (7.5) 


0 
-cos #+-cos x, 


dW Wo af sin o sin oe —— 
Hence ; _ , (7.6) 
dd 7 |COSU—COS x, COS ?—+cos Xe 
and using relation (4.2), 
dU Wo ” l ] (7 ) 
= = — . 4.d 
dg Br |cosd—cosy, cos#+cos xs 





Integrating (7.7), 


U Ue : l - tan(x,/2)—tan(3/2) 
br | sin yx, tan(y,/2)-+ tan(#/2) 


a. low COt(x »/2)—tan(#/2 \+ 4C, (7.8) 


» 
siny, ~ cot(x,/2)+ tan(3/2) 


the principal values of the logarithms being taken and C being a constant. 
Above the wing, 0 < @ < =z, and the conditions to be satisfied by wu are: 


0 6 a, u = constant = 4%,,,; 
T—Xo 6 TT, u constant Uom> for o = 0. (7.9) 


1 0 T—Xo, U 0. 
Therefore. for 0 6 ar/2, 
Ua 


hai 


Ll f ; sinh o 
tan : —yp— 
sin yx; | tan(x,/2)(cos 6+ cosh o)—sin 6 





_ if sinho ms 
tan(x,/2)(cos 8+ cosh a) +sin 6 | 


/ —sinho 
tan - —I— 
sin xo | cot(x./2)(cos 8+ cosh o)—sin 8 


a sinh o rw, 
cot( x./2)(cos 8+-cosh o)+sin 6 _ 


C0 < w/2. 


where y 0 for 0 0 XD and v = az for Ye = 
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For ml2<O0<7, 


/ ] —si f 2 
ee Dax . Die P sinh o tan(x, 2) = 
Br | sin x, | tan(x,/2)sin @— (cosh o—cos @) 


—_tan-! —sinh o tan(x,/2) \) 
, tan(y,/2)sin 6+ (cosh o—cos @) | 





— —sinh o cot(x,/2) 
— € - _ $$ J — pp — 
sin x,| cot(x./2)sin @— (cosh o—cos 6) 
heal - sini e OOM x,/2) i |! LC, 
cot(x./2)sin @+- (cosh a—cos @) 
where v = z for 7— xy, < 0 < a, and v = 0 for 7/2 < 0 < m—y,. Hence 
U x U X U xX 
C= , u ——., t. = —_——_-, 7.10 
Bsin x, im Bsinx, Sep, \ 
Therefore, on the upper surface of the wing, for 6 = 0 
2U é 2 
.= —— mites sin-! a *) + 
Bar |sinx, {tan?(y,/2)-+tanh2(o/2)}4 
/2 
a tan(x,/2) . 
sin x5 {tan*(x,/2)4 coth? 2(0/2)}+ 
2U Bs ] ‘ ] 
a Seer ae __= rexel (7.11) 
Br |sin x, (1+#)# " sin x, (1+¢?)? 
where 
t, = cot(x,/2)tanh(o/2), t, = cot(x,/2)coth(c/2), 
and for 6 = 7 
2U ] ] : l 
u = —* sin~! — + ——- sin-! —_____ |, (7.12) 
Br |sin x, (1+?) © sin x5 (1+-43)? 


where 
t, = cot(x,/2)tanh(o/2), t, = cot(x,/2)coth(o/2). 
On the under surface of the wing, w has the same numerical values, but 
the opposite signs. 
Using the relation (4.1), 


dV _ Wyif cos xy 4 COS Xp (7.13) 
d? 7 |cos#—cos x, | cos d+ cos x5 
Therefore 


tan(x, 2)—tan(d/: 2) 4 


. Ua 
V = cot x, log 
7 tan(y,/2)+tan(#/2) © 


, >) — a/9 
+ cot Xe log OMX: -) tan(v 2) 


ae D, (7.14) 
cot(xX> \+ 


)+tan(3/: 2) 
the principal values of the logarithms being taken, D being a constant. 





The 


There 


Hence 


and { 


On 


the o 


8. 
Macl 

It 
elem 


wher 


insid 


outs 


F 


insid 


outs 


ence 


7.10 


7.12 


. but 











YAWED DELTA WING AT INCIDENCE AT SUPERSONIC SPEEDS 365 


The conditions to be satisfied by v above the wing are: 


o<t< », v = constant = im? | 
T—Xo << O< 7, v = constant = a for o = 0. (7.15) 


M1 <9<7—- xX, v= 0. 
Therefore 
D Uxcot x, v, -Uaxcot x;, +Ucaecotx,. (7.16) 


t 2m 


Hence, on the upper surface of the wing, for 6 = 0, 





‘ 2U0 Oo y, sin 1 2U0 oot x2 sin=! ; 
T / t?)? 1 (1 -1,7)*? 
and for @ 7 
v 2U 0 ot x, sin . . rt —cot x2 sin ewe G 
7 : (1+ t;°)? = " (1 +t) 


On the under surface of the wing, v has the same numerical values, but 
the opposite signs. 

8. The lift of the yawed delta wing when both leading edges lie outside the 
Mach cone 

It has been shown, in § 5, that the pressure variation on a surface 
element of the wing is given by 

Ap deli pl lu, 

where p is the density of the undisturbed stream. 

Therefore, on the upper surface of the wing, for 6 = 0, 








2pU2 l 1 ] : 1 
Ap f : be ; sin 1 my + = sin-! ped (8.1) 
Br |sinx, (1+)? sinxe (1+7#,°)* 
inside the Mach cone, and 
pU« cae 
Ap a 8.2 
f Bsin x, ( ) 
outside the Mach cone. 
For ‘ed TT. 
2pU? l l ] , ] 
Ap f : ~ sin-! saat sin | (8.3) 
Bo |sinx, (1+2,°)? © sin x5 (1+-45)?. 
inside the Mach cone. and 
Ap pl be (8.4) 
Bsin x» 


outside the Mach cone 


The pressure variation on the under surface of the wing has the opposite 
sign. Therefore the pressure difference is 2\|Ap}. 
lo find the total lift, the pressure difference is integrated over the 


surface of the isosceles delta wing. 
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It can be shown that the total lift is 
2pU*ac? 


(B2—tan2A)# 


where A = (w,—w,)/2 is the angle between the chord ¢ and the z-axis 


tan(w,+w,)/2, (8.5 


Hence the lift coefficient, based on area, is given by 
, total lift 4x 


7 1 Uc? tan(w,—A) a (62—tan2)! (8.6 


C;, = 
4x sin p cosa (8.7 
= — —-— ‘ 9. / 
| cos(u+-A)cos(p—A) |* 
It is seen from (8.7) that the lift coefficient C, depends on the angle of 
yaw A, but is independent of the angle of the wing. When w, = w,, A = 0 
and the total lift is 
2pU2ac? 20U2ac? 
anthers tanw, = a a 
2 B? 608 x, 
and the lift coefficient is 40/8. 
These are the results obtained for the symmetric case. 


(8.8 


V. The yawed delta wing when one leading edge lies inside, and 
one leading edge lies outside, the Mach cone of the vertex 


9. Solution for the region outside the Mach cone 


As in §6, the velocity components u, v, w for the region outside the 


Mach cone satisfy the equation (3.3). 
On the leading edge, outside the Mach cone, 
sec y = sec x, = Btana). (9.1 


It can be shown that, for the regions outside and on the Mach cone, 


0<60<x,-—x, u= Ua/(Bsinx,), v = —Uacot x, 
—(xi—-x) < 0< 9, u = —Ua/(Bsinx,), v = Uacot x, (9.2 
X1—-x < |0| <2, w= 0, vy = 0. 


10. Solution for the region inside the Mach cone 

Equations (3.1), (4.1), (4.2) are used for the solution inside the Mac! 
cone. 

Inside the Mach cone, 0<o<@ and on the Mach cone, o =? 


-7 <.0< 7. On the leading edge of the wing, inside the Mach cone, 
secho = secho, = Btanw, (#@ = +7). 
On the wing, gd = 0, 0x6 < , 
0 7, G®%S0< 0. 


On the #-plane the Mach cone is represented by the 0-axis from @ = —7 
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to 6 -, and the upper surface of the wing inside the Mach cone by 


the two semi-infinite lines @ 0(0<a0<0),0=7(0,<0< 0). 








Fic. 6. AOB is a Mach plane of the leading edge OA of the wing. From the 
geometry of the figure, as in Fig. 3, angle ACB = x,. Angle ACD = y,— x. 











A4c F p00 
Eo 
-Tt -x, 0 xX; 
c D 6 
Fic. 7. #-plane. The points D, EZ in Fig. 7 correspond to points D, E in Fig. 8. 
- . ] 
The transformation C¢=d+ (10.1) 


cos U— COS xy 
transforms the semi-infinite region bounded by the lines @ = 0, 6= 7 
(¢ > 0);¢ = 0(0 <6 <7) in the #-plane, into the region » > 0 in the 
{-plane. The boundary is mapped on the €-axis. The values of d, k are 
chosen so that the points corresponding to (7, 0), (7,0) in the #-plane are 
the points (—k, 0), (+4, 0) in the f-plane. These values of d, k are given by: 


2d a" —_— (10.2) 
l+cosx, cosho,+cos x; 


l+cosxy, cosho,+cos x, 
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The conditions to be satisfied by w in the ¢-plane are: 
—0o <§< —k, w = 0, ) 
k<&€< +0, w= w= —Ua | 


The further transformation 


(y = 0). (10.4 


¢ = —kcosC’ (10.5 
transforms the region 7 > 0 in the ¢-plane into the semi-infinite rectangle, 
bounded by &’ = 0, yn’ = 0, &’ = zw (n’ > 0) in the @’-plane. 











con! t” 
0 a 
D . . 


Fig. 8. {’-plane. The points D, £ in Fig. 8 correspond to points D, E in Fig. 7. 


The conditions to be satisfied by w in the ¢’-plane are: 


, 


£ 0, w= 0: | 
a | (y’ > 0). (10.6 
f = #, w= W, | 

The region o > 0 (0 > 6 > —z) in the 3-plane maps into the region 


n < 0 in the ¢-plane, which is mapped into the semi-infinite rectangle 
bounded by €’ = a, 7’ 0, €° = 2n (n’ > 0). 

Therefore the region inside the Mach cone is: mapped into the basic 
semi-infinite rectangle bounded by ¢’ = 0, yn’ = 0, &’ = 2m (n’ > 0) in 
the ¢’-plane. Since { has a period 27, therefore this pattern is repeated 
throughout the region 7’ > 0 in the ¢’-plane; and dW/dZ’ is periodic in 
the ¢’-plane, with period 27. The only singularities of W or dW/d¢’ are 
at the points corresponding to the leading edge inside the Mach cone, that 
is, at the points congruent to the point (7,0) in the Z’-plane. Also the 
value of uw given by (4.2) must be finite at all points on the wing, except 
at the singularity. 

These conditions and the conditions (10.6) are satisfied by 


= Moran?’ (10.7 
dl 7 2 
Therefore 
dW wyi (cosh o,+-cos y,)! sin # cos(F/2) (10.8) 
d} —- x | cos(x,/2)(cos #—cos x,)(cos + cosh a,)! | 
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Using relation (4.2), 
dl Wet 


dd Ba 


(cosh o,+ cos x,)? cos(#/2) 





| (10.9) 


cos(x,/2)(cos #—cos x,)(cos 3 cosh o,)? 


Integrating (10.9), 


U a i 


Wt (cosh o, + COS x,)* tan r 
b7 


2? cos(x,/2)cosh?(o,/2) 


Dt of , 191 ‘vat tenons 4 € 
} l am sin(x, 2)+-(cosh o,+-cos x,)* tan r +€, (10.10) 
sin x, 2! sin(x,/2)— (cosh o,+- cos x,)* tant 


where sint = sin(7,+-7i7,) = sech(o,/2)sin(#/2), and C is a constant, the 
principal values of the logarithms being taken. 
The conditions to be satisfied by u, above the wing, are: 


0 ad Xi u constant = U,,); | a (10 11) 
x1 4G ~ WT. “ul 0 


When o = 0, 7, = 0, sinr, = (sin @/2)/(cosho,/2). Therefore when o = 0 


and x1 G TT, 


W (4-7 , 
ini ee = +e. 
Hence CU Wo/(B sin x1) Ua/(B sin x;). (10.12) 
When o = 0 and 0 < @ < y,, u = C, and therefore 
Uym = Ua/(Bsin x,). (10.13) 
On the upper surface of the wing, for 6 = 0, 
7) 0, sinh 7, = sinh(o/2)sech(o,/2), 
for ¢ 7, 
Fy 77/2, cosh 7, = cosh(a/2)sech(a,/2) (¢ > ag). 


Therefore on the upper surface of the wing we have, for 0 = 0, 


ul H l (,. of H 
u ; : — ——({Ztan : —_|—7 
Ir | COS(x,/2)cosh*(a,/2) sin x, | sin(x,/2) 
(10.14) 
siliiiges H (cosh o,-++ cos xy) ete/2) 
(cosh o,+-coshoa)?! 
and for @ TT, 
| H l fo. .f 2 
u — — 2 tan - == J— WO Te 
i | COS(y,/2)cosh*(o,/2) sin Xi\ bart 
(10.15) 


(cosh o,-++-cos yx,)* cosh(a/2) 


where H’ 


(cosh e—cosh o,)! 


[It can be verified, from (10.13), (10.14), (10.15), that w is continuous 


5092.7 


Bb 
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across the surface of the wing, and becomes infinite for o = a, that is on 
the leading edge of the wing, inside the Mach cone. 

On the under surface of the wing, wu has the same numerical values, but 
the opposite sign. 

Using relation (4.1), 


dV - | (cosh o,+-cos x,)! cos(#/2)cos 3 i 
cos( 


—— (10.16) 


X1/2)(cos #—cos x,)(cos #-+-cosh o,)! 


Therefore 


tan7+ 


y — % Sa Xi)' cosh o, 


a | 2#cos(y,/2)cosh?(c,/2) 


94 


. 2) 4 bts 
4-cot x, log sin(x, 2) (cosh o,+-Cos x,)' tant 


si "|+D, (10.17 
sin(x,/2)— (cosh o,-+cos y,)? tant 
where D is a constant; the principal values of the logarithms being taken. 
The conditions to be satisfied by v above the wing are: 
0<0<,,, >= constant = v : 
> X1 const t Vim\ foro = 0. (10.18) 
X1<9<7, v=0 
Hence 
D = w, cot a= —Uacot x;; = —Uacot x, (10.19 


Vim 


and therefore, on the upper surface of the wing, we have for @ = 0, 


v= =| ae wcbtncraie: — + cot x3 tan H )—al|, 


m |cos(x,/2)cosh?(c,/2) sin( x; 2) 
(10.20 
and for 6 = 7z, 
i Ua H’ cosh 02 + cot x; row ae ) nl]. 
z cos(x, 2)cosh?(a,/2) | sin(x;/2) 
(10.21 


On the under surface of the wing, v has the same numerical values but the 


opposite signs. 


The values of wu, v on the Mach cone (10.13), (10.19) agree with those 


found for the region outside the Mach cone (9.2). 


11. The lift of the yawed delta wing when one leading edge lies inside, ani 


one leading edge lies outside, the Mach cone 


On the upper surface of the wing the pressure difference on a surface 


element is given by: 

— a... —— —- | fotana =. )-al] 
Bx |cos(yx,/2)cosh?(a,/2) sin x, | sin(y,/2) j 

(11.1 
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inside the Mach cone, and by 


Ap = —-£— (11.2) 


outside the Mach cone. 
For 6 T, 
[72 , , 
Ap — : — Sa [2ten| 5) -7}]- 


Br | cos(x,/2)cosh*(o,/2) sin x, | sin(x,/2) 


(11.3) 
The total lift is found by integrating 2!Ap| over the surface of the 
isosceles delta wing. 
It can be shown that the total lift is 
2pU2axc* cos A /1—cot w, tanA 


} 
—s oe —Ad), 11.4 
B* cos(x,/2) B+-tana nite \ 


where c is the maximum chord of the wing and A = (w,—,)/2. 


Therefore, the lift coefficient, based on area, is 


0 total lift _ 4acosdA /(1—cotw, tana\} (11.5) 
“~~ 4pU%c? tan(w,—A) B+ cos(x,/2) B+tana : 
: sin y 4 . 

2?y cosA sin - v4 a 11.6 

x CO H sare r ore ( ) 


When w, = w, = p, A = 0, and yx, = 0, and (11.5) gives C, = 40/8. When 
w, = p, x, = 0, then 
a: BtanaA 


I 
B} 8+tanaA 


l 
} 4a cos A[tan( —A)tan p]*. (11.7) 
(11.7) gives the lift coefficient for the delta wing with one leading edge 
within the Mach cone, and one leading edge on the Mach cone. The result 
agrees with the value given in (5.8). 
When Ws fh, Wy 2A B, then 
Cr, = 4a/(B?—tan?A)!. (11.8) 
(11.8) gives the lift coefficient for the delta wing with one leading edge 
outside the Mach cone, and one leading edge on the Mach cone. The 


result agrees with the value given in (8.6). 

VI. Some numerical results for the isosceles delta wing 
12. (i) y = 45°, uw = 60°,0 <A < 45°. (Fig. 9.) 
In 0 A iS. 


D3 a 


—— cos A{ cos 2A—{cos(60° +- 2A)cos(60° — 2A)} }t ; (12.1) 
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whilst in 15° < A < 45°, 

2% 3! cosA 
[sin(75°—A)cos(60°—))]}** 


Formula (12.2) is valid up to A = 30°. 


— p=60" 


C/o == 


58 





G/ 
4 
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50 —-> 


N 














34 























Figs. 9-12. 


2=20°,.0<A< 45°. 


(ii) y 45°, pw 
In 0: A: 15°, 


(Fig. 10.) 


C,/« = 2 cos A/|cos(30°+-A)cos(30 —d)}*; 
and in 15° <A < 45°, 
C./« = 2! cos A/[sin(75°-+-A)eos(30°—A)]!. 


(ili) y = 35°, pw = 60°,0 <A < 35°. 


m6 = A= 35°. 


(Fig. 11.) 


27 cosA 


Cr/a = = 
. E’ cos 35 


* sin(25°—A)sin(25° 


| 





[cos 2A+- 4 cos 70°— 2{sin(95° —A)sin(95° +4 





(12.2) 
(12.3 
(12.4) 
LA) > 
LA)}E] (12.5) 
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whilst in 25 A = 35", 
C,/a = 2? 3? cosA] ——— _— (12.6) 
‘ sin(85° —A)cos(60° —A) 
Formula (12.6) is valid up to A = 30°. 
(iv) y = 70°, p = 30°,0 <A < 70°. (Fig. 12.) 
For 0 A < 40°, 
C,/o 2 cos A/[ cos(30° +-A)cos(30 —d)|*; (12.7) 
whilst for 40 A= 7, 
sin 70 | 
C,/0 2? cosa : a (12.8) 


sin(80° —A)cos(30° —A) 
Formula (12.8) is valid up to A = 60°. 
For A 30° in (i), A 30° in (iii), and A > 60° in (iv), the trailing edge 
of the wing lies inside the Mach cones of the base vertices. For A > 45° 
in (i) and (ii), A 35° in (ili), and A > 70° in (iv), the edge of the wing, 


for which w w,, becomes a trailing edge. 


Conclusion 

The cases when a trailing edge lies inside the Mach cone of any point 
of it require further investigation. The cases of the flat delta wing when 
one of the edges through the vertex is a trailing edge, the leading edge 
lying either inside or outside the Mach cone of the vertex, will be con- 
sidered in another paper. 

Acknowledgement is due to Professor W. G. Bickley for suggesting the 
subject of this paper, for his continued interest and help, and for many 
valuable suggestions. 
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SUMMARY 

The results developed in Part I (ref. 1) are applied to prove that the drag and 
side forces on a symmetrical wing at zero incidence are unaltered in magnitud 
when the direction of the flow is reversed, and that, for wings of restricted planform 
such that the flows over the subsonic edges are independent, the slope of the curve 
of lift against incidence is unaltered when the direction of the flow is reversed, 
provided that the Kutta—Joukowski condition is applied at the trailing edge fo 
both stream directions. These results are stated as theorems in § 2, and the possi- 
bility of obtaining more general theorems is discussed. 

It is also shown that the lift coefficient for a wing, having no subsonic edges and 
a straight trailing edge, depends only upon the Mach number, the angle of sweey 
of the trailing edge, and the average incidence over the upper and lower surfaces 
of the wing, as in (98). 


1. Introduction 

THE general linearized theory of supersonic flow past thin wings developed 
in Part I (ref. 1) is applied to prove some general results for the aero- 
dynamic forces. The notation used is the same as that of Part I and is 
not re-defined here. The numbering of the equations in this paper follows 
on consecutively from Part I. A list of the commoner symbols is given 
in an Appendix. 


2. The Flow-reversal theorems 
The two following theorems are proved in §§ 3 and 4. 


THEOREM I. Jf a thin wing, symmetrical about its chord plane, experiences 
no force normal to its chord plane when placed in a supersonic stream, then 
a reversal of the direction of the stream leaves the total force on the wing 
unaltered; the drag component, parallel to the direction of the stream, 1 
unaltered in magnitude and is reversed in direction, and the lateral com- 
ponent, normal to the stream and in the chord plane, is unaltered in magnitud: 
and direction. 

THEOREM II. Jf a thin nearly plane wing in a supersonic stream has @ 
planform such that the flows at the subsonic edges (if any) are independent, 
and if the velocity at the trailing edge satisfies the Kutta-Joukowski condition, 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 3 (1949)] 
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then a reversal of the direction of the stream leaves the slope of the curve of lift 
against incidence unaltered, provided that the Kutta—Joukowski condition is 
satisfied at the trailing edge for the new flow. 

The existence of general theorems of this kind was conjectured by 
Professor 8S. Goldstein in 1945, and the author is indebted to him for 
encouragement to prove them. 

That part of Theorem I which deals with the drag force has been stated 
before by W. D. Hayes (ref. 2) and T. von Karman (ref. 3). Hayes has 
stated a more general theorem for drag which applies to any system to 
which the linearized theory is applicable, and which includes the drag 
part of Theorem I as a special case. In general, however, the geometrical 
configurations of his systems are altered on reversing the flow. Hayes also 
stated a restricted form of Theorem II, applying only to wings having no 
subsonic edges (‘simple’ planforms in his terminology), and mentioned the 
probability of an extension to other planforms. 

Theorem II is very restricted in its scope compared with Theorem I 
and it is natural to ask whether it can be extended in any respect. This 
question does not appear to be easy to answer completely, although a 
partial answer is given here. The theorem is equivalent to a lift-reversal 
theorem for a plane wing of infinitesimal thickness whose planform satisfies 
the stated conditions, and it is shown later, in § 6, that there is no lift- 
reversal theorem for a wing of infinitesimal thickness having an arbitrarily 
warped surface, or a surface such that the local incidence is an arbitrary 
function of a single coordinate distance, being constant on lines normal 
the direction of this coordinate axis. This does not preclude the possibility 
that there is a lift-reversal theorem for some special surfaces, however, but 
these do not include the case of a linear variation of incidence. The possi- 
bility of removing the restriction on the planform has also been considered 


by the author, but without positive results. 


3. Proof of Theorem I 
For the symmetrical case of Part I, it was shown there that the drag 
force is (equation (60)) 


D 2pU? | (Fe <*o) dady. (70) 
U }z=+0 


Cz 


. « 


Ss 
Since (@4/@z)._,» = 0 on R and T, the integration may be taken over the 


whole plane z -0 instead of being limited to S. By putting 


L_(&o gies : 
soa E), - Fo(a, B), (71) 
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the potential in the plane z = +0 is 


a B 
Ff’ sft) datas : 
aan cor = ot 


—o —« 


and the expression for the drag force becomes 


= 8) x x B 
a ad a : ”m © ” # 8") dx'dp’ 
D = 4npU2B ,B) dade| 2 +. © Sola’, B’) dae ‘ 
p | | Sol B) « x a(S 5) | | V{(a- x’)(B—B’)\" ( 3 


If we assume that the wing is of finite extent, then 
fol 00, +00) = 0, (74 
and by making the further temporary assumption, for simplicity, that 
So(«, 8) is continuous everywhere, then the differentiations in equation 
(73) can be carried out to obtain 
D = 4mpU*B | | fo(a,B)dadp | | oat Sop) —, a 
J da J a Va — a \(B B yj 
where f,,, and fg denote the derivatives of f,(«, 8) with respect to a and £ 
respectively. 


—O —o —o —o 


If the flow is now reversed so that the stream velocity is directed along 
the negative x-axis, the new force, D’, in the direction of x increasing is 


| * fola’.B’) da'aB’ og 
- 4rpU2B ) dad w 
D mpU?} { { Sol, B) ¢ via (<- +S) {| v{(x’—«)(B’—B)} 


ah <sctih x 


Integration by parts gives 


’ r 8’) da’ dp’ mn 
)y — —4noI/2B ) dad j Sol x sf (77) 
. _ { | (fo,a+fo,g) dadp j Ji(o’—ayp’—By 


mae 


and by inverting the order of the double aiasiaiik we obtain finally 








x Nn a’ p’ 
, ra ew ian r (fa.+ g) dadB . 
D! = —4rpU2B , dad [ 0,07 tho 73) 
. | J Solo’. B') do‘ap | Ja aye By 


Comparison of equations (75) and (78) shows that 
DY’ = —D. (79) 
The assumption on the continuity of f, can be removed, and dis- 
continuities allowed on a finite number of arcs in the «, B-plane without 
altering this result, as may be seen by replacing the Riemann integrals 
by Stieltjes integrals. 
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Similarly, for the lateral force, Y, in the direction of the positive y-axis, 


y 2pU? | [ (: Po ¢ re) dady 
vs 


Cz oy y 0 


1 B 
; . . C ‘ . a2 ° # , 1ax'dp’ 
apl/? 8) dadB + : Sola B’) « - - é 
dnpl | | f(x, B) da | = a) | a aNE_PD (80) 


—n —a% 


If the flow is reversed, the new lateral force, Y’, is 


A / é 4 [ . fo(a’, B’) da’dp’ (81) 


tpl f(a, B) dadB| — } ; ; . 
™ i : | } Solo PL oa ap) ) | f(a’ —an(B’—B)} 


that | and by carrying out the same sequences of operations as were used for 
ation D and D’, we readily obtain the result 
y’ = ¥. (82) 
(75) J This completes the proof of Theorem I. 
: 
~ | 4. Proof of Theorem II 
Theorem II will be proved for the equivalent form mentioned in § 2, 
din | that is as a lift-reversal theorem for a plane wing of infinitesimal thickness. 
ale | Such a wing corresponds to the antisymmetrical case of Part I, with 
| \(a,8) = constant. For simplicity the incidence is taken to be such that 
- | Xa,8) = 1; no loss of generality arises since A depends only upon incidence 
i and Mach number, and these quantities remain unaltered in magnitude 
when the direction of flow is reversed. The lift on the wing is given by (69) 
of Part I, which contains only the value of the potential ¢, at points on 
the trailing edge of the wing for z = +0. 
(77 Consider the wing whose planform is shown in Fig. 3 (a) of Part I, and 
let the coordinates of the points N and P (which were left unspecified 
ly in Part I) be (a3, 8.) and («,,8,) respectively. For convenience, the potential 
(¢;),- 49 at the trailing edge is divided into three parts, which are considered 
(78 separately. These are: 
(i) a part 4’, depending upon the function G,, arising from the vortex 
sheet from the subsonic trailing edge 7, NV, which is non-zero only 
on the are 7, E of C; 
(79 (ii) a part ¢”, depending upon the function G,, arising from the vortex 
dis sheet from the subsonic trailing edge PT,, which is non-zero only 
out on the are FT, of C; 
rals (ili) a part 4”, consisting of terms given in (52) of Part I, which is 


non-zero only on the are N P of C. 
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On the are 7, VN, « = B,(B), and from (50), 








gt Pg 8) dé 
d’ == | ha lie PM = | , (5) as at (83) 
P Vi By (B)—2'} : VB, (8)- B—§} 
On the are NE = B,(8), and 
B,(B) ( By(p)—B 
’ ' 4 (x’ — , j G lé 
$! | a B) 1(€) ¢ ” ine 
vi B,(B)—2'} vB; 2(B)—B—§} 
pry v1 
From (40), when A Lk. 
Gi} x— A,(a)} : 2,/{A4(a)—A,(a)}. (85) 
Therefore, by putting € = a—A,(a) in (83) and (84), and integrating ¢’ 
along ET, with respect to y, we have 
| $' dy 
ET, 
-2 [ [1—By8)] Al) Nia a2joy) da 
| [1 =< Aa) a(a)| do 
By(B ) 
o(a)—A,( 


safc B,(B) | dB J \ry— Aj(a)| da. 


vl lea sig ant 

(86 
On inverting the order of integration, the integration with respect to f 
can be carried out, and since B,(8,) = B,(B,) and B,(A,(«)) = a on the 
are VT;, we have 


2B [ 6’ dy = 4 [ \{[Ag(a)—A,(a)][ B,(Ag(a))—a}}[1— 4 (a)] da 
ET, 1 
4 [ /{[4o(a)—Ay(a)][ Bo(Ao(a))—a]} do 


Q 
R3 





V{[ B.(B)— B,(B)[B—A,(B,(B))]} dB. (87) 
"4 
In a similar manner we can obtain the result 


R 


2B ( ¢’ dy = 4 [ \/{[ B.(8)— B,(B)|[ A2(B.(8))—B]} dB 
T:F be 
—4 [ V{[A2(a)—Ay(a )I 1— B,( ))}} dx. (88) 
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On the arc NG, from (48) with A E. 


Bop 
83 p” 2 / (PB—Axla N de, (89) 
J A \B(B)—a’) 
B,(p) 
ind on the are GP, the expression (52) with A = 1 may be expressed in 
the form 
I ; B:(B) ; 
(54 4 9 | }P “Ao )\ gy’ 9 | [{B—Al@)\ g 
A) | B,(B) x | N | B,(B)- x" 
Bil B;(B) 


4,/{[ B,(8)— B,(B) |[B—A,(B,(B))]}. (90) 
85) | where B,(f) B,| A,(B,(8))] when By < B < A,(ao). The curve a B,(B) 
can be constructed geometrically as illustrated in Fig. 4. This figure is 








ing ¢ 
ettered to correspond with the figures in Part I: 
| 
Sh | 
tof | All the straight lines are characteristic lines in the plane z = 0, a is any 
n the point on the are GP, and abcd is a characteristic parallelogram, b and c 
lying respectively on PL, LH. As a moves along GP, the locus of d is the 
curve a = B,(8) or B = A,(a) passing through H and G. That part of S 
lying to the right of « B,(8) is denoted by S,. The curve DJ and the 
rea S, are defined similarly, by changing left for right. Although the 
areas §,, S, are shown separated in Fig. 4, they may overlap if the aspect 
ratio of the wing is sufficiently large. 
From (89) and (90) we have 
- B B.(B) 
(S/ a ” a! / 4 
: , [{pb—A,(a 
2B | db” dy ,- [1— B,(B)| dB /(B nM N ae — 
. A . N | B,(B) = xf 
PN Bs B,(B) 
A 9 x B(B) 
. . n — (ax 
2 | ] B3(B)| dp / (B i(%)| dat 
, } AW \BAB)—aJ 
Bs B;(B) 
Lola 
RR > . 7 , 
{ /{[ B.(8)— B(8)|[B—A,(B.(B)) }} [1— B(B)] dB. (91) 
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After some algebra, the right-hand side of this expression may }p 
written as 


LF (Ain) Jaca) ee 


1 


+ 4 | V{[43(«) —A1(«) ][ B,(Ay(«))—a]} dat 
By(Bs) 
rs V{[a—B,(A,(«)) [ Ao(a)—A (a) |} da- 


~@ | \/{[Ao(a) —A (a) ][ Bo(Ao(a))—a]} da 


A o(ao) 
7 


ney | /{{ B.(B)— B,(B) |[B—A-2(B.(B))]} dB. (92 
Bs 


The results of equations (87), (88), and (92) may now be combined to give 


‘ As(a (B, 3)—a 
2B (4). .o dy | J (/(5 |B. eeett aif arrest b — 


TT; 
L 
4 | V{[Ag(a)—A (a) |[ B,(Ay(«))—a]} da 
H 
+ | \{(Bo(8)—By(B)|[B—A,(B,(B))]} 4B 
P 
rs 
4 | \{(As)—Ay()][>— By (Aya) 4 
L 
f 
+4 [ /{B(B)—B,(A)I.44(B.(B))—B} a8 
Ts 
Ty 
4 | V{[A3(«)—4 x) || Bo. o(x)) —a]} dx 
7 
' 
4 | \{{ B.(B) — B,(B)|[B—A,(B,(B))]} dp, (93) 
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where, for clarity, the limits of the integrations have been specified by the 
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nay | names of the points given in Fig. 4 rather than by the appropriate 
j oordinates. 

If the curve JD in Fig. 4 is given the equation « B,(B) or B = A,(a), 
in equivalent expression can be obtained (or written down immediately 
by using the formal symmetry between a and §) involving an integration 

‘over S,, ete. By taking the average of these two expressions, the last four 
integrals in (93) cancel with the corresponding terms in the equivalent 


expression, and we get the more symmetrical result 


(B—A,(x)) , (B,(B)—«\ 


2B | (¢1)2-+0 dy 1} \/ Bae J” Wl \B—Aj (a) 


) dadp-+ 


l T |{o -B(B)\ (A,(«)—B 
N \A,(a)—BJ A \a—B,(B) 





) dodB + 


« 


2 | V{[As(a)—A (0) ][ By(Ay(a))—a}} do 


“™ 2 | {[B.(8)—B,(8)|[B—A,(B,(8))]} 484 


2 /{[Ao(a«) —A4(«) ][a— B,(Ag(a))]} dat 





| 2 | {LB.(8)— By(B)][4(Bx(8))—B]} a8. (94) 


This expression now possesses complete formal symmetry, not only with 
respect to a and £, but also with respect to the fore and aft on the wing, 
so that a reversal of the stream direction merely reverses the sign of the 


lift force without altering its magnitude. This proves Theorem II. 


5. The lift of a wing with a straight supersonic trailing edge and 
no subsonic edges 
Before considering the possibility of a more general theorem for lift than 
Theorem II, it is advantageous to obtain an interesting result for the lift 
loree on wings which have no subsonic edges, and whose trailing edge is 
straight. 
(93 Since the edges of S are supersonic, the flows over the two surfaces are 
independent. The potential due to an element of area dS at (x,y, +0) on 


Sat a point (x,,y¥,, +0) on the trailing edge and inside the downstream 


y the 
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characteristic cone with its vertex at (x,y, 0) is, from (7) and (13) of Part 
l, dS 
dd — r . 5) 2)° (95 
m7 {(%,—2)?— By, —y)*} 

The lift, dL, from this element is obtained by integrating dd with 
respect to y, along the relevant portion of the trailing edge. If 0 is the 
angle of sweep of the trailing edge, it is easy to show that 
l,cosQ dS 


iL, = pU* . 
: P / (M* cos? — 1) 


(96 


A similar result is obtained for the lower surface. Thus the total lift, L. js 


pU*cosO a : 
an (l,—1,) dS, (97 
/(M?cos*?O—1) J , 
s 
and the lift coefficient with respect to }pU? and the wing area is 
4cos® —— 
se 3(/,—I,), (98 
/ (M? cos? — 1) 
where the bar denotes the average value with respect to the area. The 
quantity 4(/,—/,) is the average value of the local incidences on the upper 
and lower surfaces, if the incidence is defined as a right-hand rotation 
about the y-axis. 

Special cases of this result have been noticed before, usually in con- 
nexion with plane triangular wings. The special case of this formula (98 
for a plane wing is a consequence of Theorem II and the well-known result 
for the lift coefficient of a ‘two-dimensional’ swept wing, due to Busemann. 


6. The possibility of a more general theorem than Theorem II 
By using the result given in (96) it can be shown that there is no general 
lift-reversal theorem for a wing of infinitesimal thickness having an arbi- 
trarily warped surface. For if there is such a theorem it must be true for 
special cases, and, in particular, it must be true for a wing having the plan- 
form shown in Fig. 5. If the local incidence is varied over an area d8 at 
the point (x,y) shown, which is such that the characteristic lines through 
it cut only the straight leading and trailing edges, then equation (9%) 
shows that the increments of lift due to this variation are unequal for 
the two directions of flow, since the angles of sweep of the two supersonic 
edges are unequal. Thus the theorem cannot be true both before and after 
the incidence variation, and so there can be no general theorem of this kind. 
In a similar way, it can be argued that there is no general lift-reversal 
theorem for variations of incidence such that the incidence is constant 
along any set of parallel straight lines and varies arbitrarily in a direction 
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normal to them. In particular, the incidence cannot be an arbitrary 


function of x alone, or of y alone. 








¥ 
dS K(x, y) 
y + / \ 
he 
Fic. 5. 


For any planform with no subsonic edges, we can write down 
A Bo B. (8) B . 

2B [ (?,).- +0 dy | ] 4(B) ]4p | 3 ih , | A(a' x’, B’) dp’ 

| : ( (B)—a' \(B—B’) 

| (99) 


T2T1 0 By(B) Aj(a’) 
which, after some transformation, 


3—A,(a)) | ( B,(B)—«) 
B,(B)—a) © J \B—A,(a)} 


{ (, eet “f + /{AHe |) dadB. (100) 


A;(a) 





Se) 


Jala, Au(o))4 


js) 


The corresponding quantity when the flow is reversed is 


~t{ (P—A(o)) Faagae| B = 
4 (, | B,(B) x} AJ \B—A,(a)} ( 2(B), 8) 





B) 

F x —A,(x) 

| » ie (8 (9) OX aa dodB, (101) 
N |B—A,(a)} | a’—a ); 

and it can be verified that, if A is any linear function of « and f, the 

magnitudes of these two quantities are unequal in general. Thus 

appears probable that Theorem II cannot be generalized in respect of 


the form of the surface for a general planform. 

If we restrict the planform to be symmetrical about the z-axis, then 
| we have at once the almost trivial result that a lift-reversal theorem holds 
for an incidence variation with y alone and odd in y, since the lift vanishes 
identically, by symmetry. There appears to be no corresponding result 
when the incidence is a function of x. 

The problem of removing the restriction on the planform has not been 
solved. The author’s personal opinion is that some more general theorem 
exists for restricted classes of planforms, but that Theorem II is all that 
is true for general planforms. 
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Since this paper was written, the theorems of § 2 have been stated by 
Hayes (ref. 4). The author has not seen ref. 4, but he understands that 
Mr. Hayes has found a counter-example to show that Theorem II is the 


most general theorem for general planforms. 
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APPENDIX: LIST OF SYMBOLS 
A,(a), Ag(ax) functions defining C, equation (21). 
B,(B), B,(B) functions defining C, equation (22). 
B ,/(M?2—1). 
BY, L rectangular components of the aerodynamic force on the 
wing. 
G,, G, functions depending on the strength of the vortex sheet, 
equations (40) and (51). 
Mach number. 
velocity of the undisturbed stream. 


1 (ad 
“oa. 


direction cosines of outward normal to upper wing sur- 
face (z +-(), 


+-0° 


direction cosines of outward normal to lower wing sur- 
face (27 = —0). 
rectangular components of the disturbance velocity, 
U grad d. 
rectangular Cartesian coordinates. 
projection of the wing planform boundary on z = 0. 
regions of the plane z = 0 (ef. § 2, Part I). 
x, B x+ By, characteristic coordinates in z = 0. 
(a, 9), (0,85), (a458:) \ characteristic coordinates of the points L, M, T,, T,, N, P 
(cas Be)s (a%3»B3)s (4s Bs)! on C respectively (Fig. 1, Part I, or Fig. 4). 
Xa; B) value of f(a,8) on S. 
je x, B), v( x, B) values of f(a, B) on R. 
s(x; B), v(x, B) values of f(a, B) on T. 
d potential function such that (u,v, w) U grad ¢. 
dos parts of ¢ respectively even and odd functions of z. 
p density of the undisturbed stream. 
Other symbols which do not occur frequently in the text are defined wherever 
they are used. 
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